{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Content under Creative Commons Attribution license CC-BY 4.0, code under BSD 3-Clause License © 2018 D. Koehn, notebook style sheet by L.A. Barba, N.C. Clementi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href=\"https://fonts.googleapis.com/css?family=Merriweather:300,300i,400,400i,700,700i,900,900i\" rel='stylesheet' >\n",
       "<link href=\"https://fonts.googleapis.com/css?family=Source+Sans+Pro:300,300i,400,400i,700,700i\" rel='stylesheet' >\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' >\n",
       "<style>\n",
       "\n",
       "@font-face {\n",
       "    font-family: \"Computer Modern\";\n",
       "    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "}\n",
       "\n",
       "\n",
       "#notebook_panel { /* main background */\n",
       "    background: rgb(245,245,245);\n",
       "}\n",
       "\n",
       "div.cell { /* set cell width */\n",
       "    width: 800px;\n",
       "}\n",
       "\n",
       "div #notebook { /* centre the content */\n",
       "    background: #fff; /* white background for content */\n",
       "    width: 1000px;\n",
       "    margin: auto;\n",
       "    padding-left: 0em;\n",
       "}\n",
       "\n",
       "#notebook li { /* More space between bullet points */\n",
       "margin-top:0.5em;\n",
       "}\n",
       "\n",
       "/* draw border around running cells */\n",
       "div.cell.border-box-sizing.code_cell.running { \n",
       "    border: 1px solid #111;\n",
       "}\n",
       "\n",
       "/* Put a solid color box around each cell and its output, visually linking them*/\n",
       "div.cell.code_cell {\n",
       "    background-color: rgb(256,256,256); \n",
       "    border-radius: 0px; \n",
       "    padding: 0.5em;\n",
       "    margin-left:1em;\n",
       "    margin-top: 1em;\n",
       "}\n",
       "\n",
       "\n",
       "div.text_cell_render{\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    line-height: 140%;\n",
       "    font-size: 110%;\n",
       "    width:680px;\n",
       "    margin-left:auto;\n",
       "    margin-right:auto;\n",
       "}\n",
       "\n",
       "/* Formatting for header cells */\n",
       ".text_cell_render h1 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-style:regular;\n",
       "    font-weight: bold;    \n",
       "    font-size: 250%;\n",
       "    line-height: 100%;\n",
       "    color: #004065;\n",
       "    margin-bottom: 1em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       ".text_cell_render h2 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: bold; \n",
       "    font-size: 180%;\n",
       "    line-height: 100%;\n",
       "    color: #0096d6;\n",
       "    margin-bottom: 0.5em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       "\n",
       ".text_cell_render h3 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "\tfont-size: 150%;\n",
       "    margin-top:12px;\n",
       "    margin-bottom: 3px;\n",
       "    font-style: regular;\n",
       "    color: #008367;\n",
       "}\n",
       "\n",
       ".text_cell_render h4 {    /*Use this for captions*/\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: 300; \n",
       "    font-size: 100%;\n",
       "    line-height: 120%;\n",
       "    text-align: left;\n",
       "    width:500px;\n",
       "    margin-top: 1em;\n",
       "    margin-bottom: 2em;\n",
       "    margin-left: 80pt;\n",
       "    font-style: regular;\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {  /*Use this for small titles*/\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    font-weight: regular;\n",
       "    font-size: 130%;\n",
       "    color: #e31937;\n",
       "    font-style: italic;\n",
       "    margin-bottom: .5em;\n",
       "    margin-top: 1em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h6 { /*use this for copyright note*/\n",
       "    font-family: 'Source Code Pro', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 9pt;\n",
       "    line-height: 100%;\n",
       "    color: grey;\n",
       "    margin-bottom: 1px;\n",
       "    margin-top: 1px;\n",
       "}\n",
       "\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "\t\n",
       "    \n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"], \n",
       "                           equationNumbers: { autoNumber: \"AMS\", useLabelIds: true}\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Execute this cell to load the notebook's style sheet, then ignore it\n",
    "from IPython.core.display import HTML\n",
    "css_file = '../style/custom.css'\n",
    "HTML(open(css_file, \"r\").read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Equidistant Cartesian grid\n",
    "\n",
    "In [lecture 1](https://github.com/daniel-koehn/Theory-of-seismic-waves-II/tree/master/01_Analytical_solutions) we reviewed the forward problems to describe wave propagation in isotropic elastic and acoustic media. Furthermore, analytical solutions for homogeneous 3D, 2D and 1D acoustic media are derived to verify the accuarcy of the finite difference codes we develop later in the course.\n",
    "\n",
    "Before introducing the basics of finite-difference modelling, I want to give an overview of some concepts to generate discrete subsurface models."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## From the continous to discrete world  \n",
    "\n",
    "The partial differential equations, which we developed before describe wave propagation in a continous medium. However, if we want to solve these equations in a computer, we have to discretize the continuum, due to the limitation that a computer can only compute discrete problems. The question is how fine do we have to discretize the continuum in order to get accurately seismic modelling results?\n",
    "\n",
    "Do we have to discretize down to the molecular scale? For example 12 grams of the carbon isotop C-12 contains $6.022140857x10^{23}$ atoms. Let's assume that we can represent the coordinates of the atoms in each direction by a double precision floating-point number which requires 8 Byte, so in totally we need ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Memory required =  14453138056800.0  TB\n"
     ]
    }
   ],
   "source": [
    "memory_mole_c12 = 6.022140857e23 * 3 * 8 * 1e-12\n",
    "print (\"Memory required = \", memory_mole_c12, \" TB\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$1.4x10^{13}$ TB would lead to some serious memory problems when using such a fine mesh, because the [NEC-cluster](https://www.rz.uni-kiel.de/en/our-portfolio/hiperf/nesh) at the computing centre at Christian-Albrechts-University Kiel only has a total of 100 TB of RAM. We have not even considered the required runtime for such a seismic modelling run.\n",
    "\n",
    "Fortunately, we do not have to discretize the universe down to the atomic or moleculare scale in order to get reasonable modelling results. In a later lecture, we will derive criteria to find the best compromise between accuracy of the modelling results and level of model discretization."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discretization of a 2D model\n",
    "\n",
    "Let's start with the 2D acoustic wave equation:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2 P(x,z,t)}{\\partial t^2} - Vp^2(x,z) \\biggl(\\frac{\\partial^2 P(x,z,t)}{\\partial x^2} + \\frac{\\partial^2 P(x,z,t)}{\\partial z^2} \\biggl)= f \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "In a first step to solve this problem using finite-differences we have to discretize the pressure wavefield $P(x,z,t)$ and $V_{p0}(x,z)$ at a specific time t on a **equidistant 2D Cartesian spatial grid**:\n",
    "\n",
    "<img src=\"images/2D-grid_cart_ac.png\" width=\"75%\">\n",
    "\n",
    "where we define the positions (x,z) of the discrete points in the model by\n",
    "\n",
    "\\begin{align}\n",
    "x &= i*dh,\\nonumber \\\\\n",
    "z &= j*dh,\\nonumber \\\\\n",
    "\\end{align}\n",
    "\n",
    "where $dh$ is the spatial grid point distance. We denote the value of the pressure wavefield and P-wave velocity at a discrete point by $P_{j,i}$ and $Vp_{j,i}$, respectively. These are 2D matrices, so we can simply store them in 2D arrays. Let's try to do this in Python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Import Libraries\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Here, I introduce a new library, which is useful \n",
    "# to define the fonts and size of a figure in a notebook\n",
    "from pylab import rcParams\n",
    "\n",
    "# Get rid of a Matplotlib deprecation warning\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to discretize a homogeneous subsurface model with x = 1000 m and z = 500 m and a spatial grid point distance dh = 50 m."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "# Define FD grid parameters\n",
    "dh = 50.0\n",
    "Xmax = 1000.0\n",
    "Zmax = 500.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define no. of grid points in x- and z-direction and grid point coordinates\n",
    "def coord_def(Xmax,Zmax,dh):\n",
    "    \n",
    "    # number of grid points in x- and z-direction\n",
    "    NX = (int)(Xmax / dh)\n",
    "    NZ = (int)(Zmax / dh)\n",
    "    \n",
    "    print(\"NX = \",NX)\n",
    "    print(\"NZ = \",NZ)\n",
    "    \n",
    "    # x- and z-coordinates of the discrete Cartesian grid points\n",
    "    x = np.arange(0, (NX+1)*dh, dh)\n",
    "    z = np.arange(0, (NZ+1)*dh, dh)\n",
    "    \n",
    "    XX, ZZ = np.meshgrid(x, z)\n",
    "    \n",
    "    return NX, NZ, x, z, XX, ZZ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NX =  20\n",
      "NZ =  10\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtMAAAFNCAYAAADCcOOfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3df5yld13f/deHhEw4EElCDCE/OEOKtVDqj8wIS8ttp4IWIgp9FHXFlog/YoU+sHdLS2wfCnjfW4V6+4NbKipQbSeQRNQKUYoxsFrrvYFZAiHcgTs/2N3EbGBjfoCAI4TP/cd1DTsze2Z29trz87Ov5+NxPfac7znn874+c+ZkP7nmumYjM5EkSZJ04h416R2QJEmSZpXDtCRJktSRw7QkSZLUkcO0JEmS1JHDtCRJktSRw7QkSZLUkcO0pFNGRDw5Iv4qIk7b4vHXRcTyTp6r8YmI90bEFZPeD9h+XyJiPiIyIk4f935JmhyHaUlTJyIORMQX22F2bfuVk62bmYcy83GZ+ciwnhsRPxgRf3ay+7au3lcH+i0en4uIt0XEwYj4XETcHBEvWPf4UkR8Zd3X7Z6IuC4ivuU4uWe02bdHxOfb9+DtETHfsY+liLiny2s3y8wXZOZvDaPWyZqmfZE0HRymJU2r72qH2bXtX056h6bE6cDdwD8EHg/8FHDdpqH33sx8HHAWsAv4BPA/I+K529R9F/DdwEvbut8I7Ae2e81AFY/MRsO/MyUdw/8wSJopEXFaRPx8RNwfEXdFxCvX/2i9PaL6vHXPX3/qxvym5z4lIv6kPcJ7A3Deutdtfu4Ptnmfi4hPRcQPRMTTgLcAz26PAj/UPvc72yPGn42IuyPidQPqXhERh9o+/kP72POBfw98X1vvo5v7z8zPZ+brMvNAZn4lM68HPgUsDHhuZuY9mfnTwFuBN2zxNX0e8O3AizLzQ5n55cx8ODPfnJlva5/z8oi4re3/roj4sXWvX2qPgL8mIu4D3gm8F7hw3RHyCyPiURFxVUTcGRF/2R4xP7etcWZELLfrD0XEhyLiie1jeyPiR9rbfysi3t8+7/6IuDoizl63Lwci4tURcUtEPBwR10bEmVv0fVpE/F9tnU9FxL/c9J7vjYg9EfG/gC8Al27alw3fi8B3DsqRVJvDtKRZ86PAC4FvBhaBl5xErXfQHH09D/g/gK3OhX0s8CbgBZl5FvD3gY9k5m3AvwD+n/bo+dpQ93ngZcDZNAPWj0fEizeVfQ7w9TRHfn86Ip6Wmf8D+I/AtW29bzxeA+3A+beBjx/nqb8LXNb2stnzgA9m5t3bvP4zNF/3rwFeDvxiRFy27vELgHOBPk3vL6A9Qt5u9wKvAl5Mc1T9QuBB4M3t66+gOSJ+CfAEmq/rFwe1DPxs+/qntc9/3abnfC/wfOApwDcAP7hFTz/a7uc3AZe1+7bZPweupDnKf3DA64f1vShpRjlMS5pW/709Qrm2/Wi7/r3AL2Xm3Zn5AM1gdcIi4snAtwA/lZmrmfmnwHu2eclXgGdExGMy83Bmbjm8ZubezPxYe+T4Fpojtf9w09Nen5lfzMyPAh+lOa3iRHt4NHA18FuZ+YnjPP1emkH07AGPPQE4vN2LM/MPMvPO9mj3nwB/BPxv657yFeC17ddy0BAM8GPAf2iPlq/SDMEvaY8Ef6ndj6dm5iOZuT8zPztgP+7IzBvanCPAL3Ds1/ZNmXlv+/3xHppheZDvBX653Z8HgZ8b8JzfzMyPt0frvzTg9Sf9vShptjlMS5pWL87Ms9dtv9GuX0hzzvCazUcLd+pC4MHM/PzxarXP+T6ao6WHI+IPIuLvbFU4Ip4VER+IiCMR8XD7uvM2Pe2+dbe/ADzuRHa+PX/3vwF/A+zkfPKLgAQeGvDYXwJPOk7eCyJiX0Q80J7OcjkbezqSmX99nH3oA7+39j9IwG3AI8AT217eB1wTEfdGxBvb/1nYvB/nR8Q1EfEXEfFZYJnuX9vN30uDjsxvd7R+WN+LkmaYw7SkWXOY5kf7a5686fHPA7119y/Yps45m0572FzrqzLzfZn57TRD5yeAteE+Bzz9HcC7gUsy8/E051XHVrU3Rx3vCRERwNtohtB/OuCI6SD/BPjwpv95WPPHwDMj4uIt8uaA3wF+HnhiezrLH7Kxp837PaiPu2lOlVn/P0lnZuZfZOaXMvP1mfl0mtNoXkhzushmP9vW/obM/Brgn7Hzr+1mh4H1PV8y4DnbvR/H+16UdApwmJY0a64DXhURF0fEOcBVmx7/CLA7Ih4dEVuex5qZB4EV4PXR/Fq45wDfNei5EfHEiPjudvBeBf6K5ogqwKeBiyPijHUvOQt4IDP/OiKeSfMbMnbq08B8bP+bI36V5nzh79rmlIq130BxUUS8FvgRmosbj5GZfwzcQHPUeCEiTo+IsyLiX0TEDwFnAHPAEeDL0fwqvu/YQR9PiIjHr1t7C7AnIvrt/n1tRLyovf2PIuLvRfN7vT9Lc9rHoF9LeBbN1/+hiLgI+LfH2Y/tXAf8RPs1Oht4TYfXb/e9KOkU4DAtaVq9Jzb+nunfa9d/g+Z0gI8CH6a5sG69nwL+Fs3Fba+nOUq8lZcCzwIeAF4L/Nctnvco4N/QnHf8AM05uq9oH3s/zcV/90XE/e3aK4CfiYjPAT9NM3Tt1G+3f/5lRHx484PtIPpjNOcB37fu6/MD6552YUT8Fc3Q+SHg7wFLmflH2+S+hOZo87XAw8CtNBfV/XFmfo7m4sHraL6uL6U58r6l9hzudwJ3tad1XAj8cvu6P2q/Nvtovv7Q/AThXTSD9G3An9CcwrHZ62kuFnwY+AOOff9PxG/QnPt9C3AzTf9fZvAQv9Xrt/telHQKiMzj/kRRkqZWNL9f+VPAozPzy5PdG82y9oj7WzKzP+l9kTQ7PDItSTolRcRjIuLy9rSWi2h+OvF7x3udJK3nMC1JOlUFzWkjD9Kc5nEbzWk5krRjnuYhSZIkdeSRaUmSJKkjh2lJkiSpo9MnvQMn47zzzsv5+flJ74YkSZKK279///2Z+bWb12d6mJ6fn2dlZWXSuyFJkqTiIuLgoHVP85AkSZI6cpiWJEmSOnKYliRJkjpymJYkSZI6cpiWJEmSOnKYliRJkjpymJYkSZI6mqphOiKeHxGfjIg7IuKqSe/Pqerqq2F+Hh71qObPq6+erfpmTE/9KhkVeqiSUaGHcWRU6KFKRoUedByZORUbcBpwJ3ApcAbwUeDp271mYWEhNVzLy5m9XiYc3Xq9Zn0W6psxPfWrZFTooUpGhR7GkVGhhyoZFXrQUcBKDphHo3ls8iLi2cDrMvMft/d/EiAzf3ar1ywuLqb/AuJwzc/DwQH/vs/c3H3s2rX7pOvv23cNq6sXjKy+GdNTv0pGhR6qZFToYRwZFXqokjHJHvp9OHBgKBFqRcT+zFzcvD5Np3lcBNy97v497doGEXFlRKxExMqRI0fGtnOnikOHBq+vrp4/lPpb1RlWfTOmp36VjAo9VMmo0MM4Mir0UCVjkj1s9fe5RmDQ4epJbMD3AG9dd/+fA//3dq/xNI/h6/c3/qhobev3Z6O+GdNTv0pGhR6qZFToYRwZFXqoklGhBx3FFqd5TNOR6XuAS9bdvxi4d0L7csraswd6vY1rvV6zPgv1zZie+lUyKvRQJaNCD+PIqNBDlYwKPWgHBk3Yk9iA04G7gKdw9ALEv7vdazwyPRrLy5lzc4cTHsl+f/gXMYy6vhnTU79KRoUeqmRU6GEcGRV6qJJRoQc1mPYLEAEi4nLgl2h+s8fbM3Pb/6/yAsTRWVpaAmDv3r0zWd+M6alfJaNCD1UyKvQwjowKPVTJqNCDtr4A8fRJ7MxWMvMPgT+c9H5IkiRJOzFN50xLkiRJM8VhWpIkSerIYVqSJEnqyGFakiRJ6shhWpIkSerIYVqSJEnqyGFakiRJ6shhWpIkSerIYVqSJEnqyGFakiRJ6shhWpIkSerIYVqSJEnqyGFakiRJ6shhWpIkSerIYVqSJEnqyGFakiRJ6iozZ3ZbWFhIDd/ycubc3OGER7Lfb+7PUn0zpqd+lYwKPVTJqNDDODIq9FAlo0IPagArOWAenfhAfDKbw/TwLS9n9nrNd8ba1usN74M56vpmTE/9KhkVeqiSUaGHcWRU6KFKRoUedNRWw3Q0j82mxcXFXFlZmfRulDI/DwcPHrs+N3cfu3btPun6+/Zdw+rqBSOrb8b01K+SUaGHKhkVehhHRoUeqmRMsod+Hw4cGEqEWhGxPzMXN697zrQ2OHRo8Prq6vlDqb9VnWHVN2N66lfJqNBDlYwKPYwjo0IPVTIm2cNWf59rBAYdrp6VzdM8hq/f3/ijorWt35+N+mZMT/0qGRV6qJJRoYdxZFTooUpGhR50FFuc5uGRaW2wZw/0ehvXer1mfRbqmzE99atkVOihSkaFHsaRUaGHKhkVetAODJqwZ2XzyPRoeGXzqZNRoYdxZFTooUpGhR7GkVGhhyoZFXpQAy9A1IlYWloCYO/evTNZ34zpqV8lo0IPVTIq9DCOjAo9VMmo0IO8AFGSJEkaOodpSZIkqSOHaUmSJKkjh2lJkiSpI4dpSZIkqSOHaUmSJKkjh2lJkiSpI4dpSZIkqSOHaUmSJKkjh2lJkiSpI4dpSZIkqSOHaUmSJKkjh2lJkiSpI4dpSZIkqSOHaUmSJKkjh2lJkiSpo5EN0xHx9oj4TETcum7t3Ii4ISJub/88p12PiHhTRNwREbdExGWj2i9JkiRpaDJzJBvwrcBlwK3r1t4IXNXevgp4Q3v7cuC9QAC7gJt2krGwsJAavuXlzLm5wwmPZL/f3J+l+mZMT/0qGRV6qJJRoYdxZFTooUpGhR7UAFZy0Mw7aHFYGzC/aZj+JPCk9vaTgE+2t38N+P5Bz9tuc5gevuXlzF6v+c5Y23q94X0wR13fjOmpXyWjQg9VMir0MI6MCj1UyajQg47aapiO5rHRiIh54PrMfEZ7/6HMPHvd4w9m5jkRcT3wc5n5Z+36jcBrMnNlu/qLi4u5srLtU3SC5ufh4MFj1+fm7mPXrt0nXX/fvmtYXb1gZPXNmJ76VTIq9FAlo0IP48io0EOVjEn20O/DgQNDiVArIvZn5uLm9Wm5ADEGrA2c8iPiyohYiYiVI0eOjHi3Tj2HDg1eX109fyj1t6ozrPpmTE/9KhkVeqiSUaGHcWRU6KFKxiR72Orvc43AoMPVw9rwNI+Z0+9v/FHR2tbvz0Z9M6anfpWMCj1UyajQwzgyKvRQJaNCDzqKLU7zGPeR6XcDV7S3rwB+f936y9rf6rELeDgzD4953wTs2QO93sa1Xq9Zn4X6ZkxP/SoZFXqoklGhh3FkVOihSkaFHrQDgybsYWzAO4HDwJeAe4AfBp4A3Ajc3v55bvvcAN4M3Al8DFjcSYZHpkfDK5tPnYwKPYwjo0IPVTIq9DCOjAo9VMmo0IMaTOICxFHzAsTRWVpaAmDv3r0zWd+M6alfJaNCD1UyKvQwjowKPVTJqNCDpv8CREmSJGnmOExLkiRJHTlMS5IkSR05TEuSJEkdOUxLkiRJHTlMS5IkSR05TEuSJEkdOUxLkiRJHTlMS5IkSR05TEuSJEkdOUxLkiRJHTlMS5IkSR05TEuSJEkdOUxLkiRJHTlMS5IkSR05TEuSJEldZebMbgsLC6nhW17OnJs7nPBI9vvN/Vmqb8b01K+SUaGHKhkVehhHRoUeqmRU6EENYCUHzKMTH4hPZnOYHr7l5cxer/nOWNt6veF9MEdd34zpqV8lo0IPVTIq9DCOjAo9VMmo0IOO2mqYjuax2bS4uJgrKyuT3o1S5ufh4MFj1+fm7mPXrt0nXX/fvmtYXb1gZPXNmJ76VTIq9FAlo0IP48io0EOVjEn20O/DgQNDiVArIvZn5uLmdc+Z1gaHDg1eX109fyj1t6ozrPpmTE/9KhkVeqiSUaGHcWRU6KFKxiR72Orvc43AoMPVs7J5msfw9fsbf1S0tvX7s1HfjOmpXyWjQg9VMir0MI6MCj1UyajQg45ii9M8PDKtDfbsgV5v41qv16zPQn0zpqd+lYwKPVTJqNDDODIq9FAlo0IP2oFBE/asbB6ZHg2vbD51Mir0MI6MCj1UyajQwzgyKvRQJaNCD2rgBYg6EUtLSwDs3bt3JuubMT31q2RU6KFKRoUexpFRoYcqGRV6kBcgSpIkSUPnMC1JkiR15DAtSZIkdeQwLUmSJHXkMC1JkiR15DAtSZIkdeQwLUmSJHXkMC1JkiR15DAtSZIkdeQwLUmSJHXkMC1JkiR15DAtSZIkdeQwLUmSJHXkMC1JkiR15DAtSZIkdeQwLUmSJHU0smE6Ii6JiA9ExG0R8fGI+Il2/dyIuCEibm//PKddj4h4U0TcERG3RMRlo9o3SZIkaSgycyQb8CTgsvb2WcD/BzwdeCNwVbt+FfCG9vblwHuBAHYBNx0vY2FhITV8y8uZc3OHEx7Jfr+5P0v1zZie+lUyKvRQJaNCD+PIqNBDlYwKPagBrOSgmXfQ4ig24PeBbwc+CTwpjw7cn2xv/xrw/eue/9XnbbU5TA/f8nJmr9d8Z6xtvd7wPpijrm/G9NSvklGhhyoZFXoYR0aFHqpkVOhBR201TEfz2GhFxDzwp8AzgEOZefa6xx7MzHMi4nrg5zLzz9r1G4HXZObKVnUXFxdzZWXLh9XB/DwcPHjs+tzcfezatfuk6+/bdw2rqxeMrL4Z01O/SkaFHqpkVOhhHBkVeqiSMcke+n04cGAoEWpFxP7MXNy8PvILECPiccDvAP8qMz+73VMHrB0z6UfElRGxEhErR44cGdZuqnXo0OD11dXzh1J/qzrDqm/G9NSvklGhhyoZFXoYR0aFHqpkTLKHrf4+1wgMOlw9rA14NPA+4F+vW/M0jynW72/8UdHa1u/PRn0zpqd+lYwKPVTJqNDDODIq9FAlo0IPOootTvMY5W/zCOBtwG2Z+QvrHno3cEV7+wqac6nX1l/W/laPXcDDmXl4VPunwfbsgV5v41qv16zPQn0zpqd+lYwKPVTJqNDDODIq9FAlo0IP2oFBE/YwNuA5NKdp3AJ8pN0uB54A3Ajc3v55bvv8AN4M3Al8DFg8XoZHpkfDK5tPnYwKPYwjo0IPVTIq9DCOjAo9VMmo0IMaTPICxFHxAsTRWVpaAmDv3r0zWd+M6alfJaNCD1UyKvQwjowKPVTJqNCDJngBoiRJklSVw7QkSZLUkcO0JEmS1JHDtCRJktSRw7QkSZLUkcO0JEmS1JHDtCRJktSRw7QkSZLUkcO0JEmS1NHp2z0YEZ89zusDOJyZf3t4uyRJkiTNhm2HaeDOzPzm7Z4QETcPcX8kSZKkmXG80zz+6Q5q7OQ5kiRJUjnbHpnOzLvW34+Ir1n/msx8YPNzJEmSpFPF8U7zACAifgz4GeCLQLbLCVw6ov2SJEmSpt6Ohmng1cDfzcz7R7kzkiRJ0izZ6a/GuxP4wih3RJIkSZo1Oz0y/ZPAn0fETcDq2mJmvmokeyVJkiTNgsw87gZ8EPgF4OXAFWvbTl47ym1hYSE1fMvLmXNzhxMeyX6/uT9L9c2YnvpVMir0UCWjQg/jyKjQQ5WMCj2oAazkoDl50OIxT4I/38nzxr05TA/f8nJmr9d8Z6xtvd7wPpijrm/G9NSvklGhhyoZFXoYR0aFHqpkVOhBR201TEfz2PYiYg9wEHgPG0/zeGDIB8pPyOLiYq6srExyF8qZn4eDB49dn5u7j127dp90/X37rmF19YKR1TdjeupXyajQQ5WMCj2MI6NCD1UyJtlDvw8HDgwlQq2I2J+Zi5vXd3rO9EvbP39y3Zq/Gq+gQ4cGr6+unj+U+lvVGVZ9M6anfpWMCj1UyajQwzgyKvRQJWOSPWz197lGYNDh6lnZPM1j+Pr9jT8qWtv6/dmob8b01K+SUaGHKhkVehhHRoUeqmRU6EFHscVpHtv+aryIuOx4w/hOnqPZsWcP9Hob13q9Zn0W6psxPfWrZFTooUpGhR7GkVGhhyoZFXrQDgyasNc24KPAOcC522w3b1djlJtHpkfDK5tPnYwKPYwjo0IPVTIq9DCOjAo9VMmo0IMadLkAMSIOAF8BYpt5/EhmPnMIc/0J8wLE0VlaWgJg7969M1nfjOmpXyWjQg9VMir0MI6MCj1UyajQgzpegJiZ8yPbI0mSJGnG7fSfE5ckSZK0icO0JEmS1JHDtCRJktTRjobpiLgxIi7ftPbro9klSZIkaTbs9Mj0U4DXRMRr160dczWjJEmSdCrZ6TD9EPBc4IkR8Z6IePwI90mSJEmaCTsdpiMzv5yZrwB+B/gzYHj/sLwkSZI0g7b9PdPrvGXtRmb+ZkR8DHjlaHZJkiRJmg07GqYz89c23d8P/NBI9kiSJEmaEf5qPEmSJKkjh2lJkiSpI4dpSZIkqSOHaUmSJKkjh2lJkiSpI4dpSZIkqaORDdMRcWZEfDAiPhoRH4+I17frT4mImyLi9oi4NiLOaNfn2vt3tI/Pj2rfJEmSpKHIzJFsQACPa28/GrgJ2AVcB+xu198C/Hh7+xXAW9rbu4Frj5exsLCQGr7l5cy5ucMJj2S/39yfpfpmTE/9KhkVeqiSUaGHcWRU6KFKRoUe1ABWctDMO2hx2BvQAz4MPAu4Hzi9XX828L729vuAZ7e3T2+fF9vVdZgevuXlzF6v+c5Y23q94X0wR13fjOmpXyWjQg9VMir0MI6MCj1UyajQg47aapiO5rHRiIjTgP3AU4E3A/8J2JeZT20fvwR4b2Y+IyJuBZ6fmfe0j90JPCsz79+q/uLiYq6srIxs/09F8/Nw8OCx63Nz97Fr1+6Trr9v3zWsrl4wsvpmTE/9KhkVeqiSUaGHcWRU6KFKxiR76PfhwIGhRKgVEfszc3Hz+kgvQMzMRzLzm4CLgWcCTxv0tPbP2Oaxr4qIKyNiJSJWjhw5MrydFQCHDg1eX109fyj1t6ozrPpmTE/9KhkVeqiSUaGHcWRU6KFKxiR72Orvc43AoMPVo9iA1wL/Fk/zmGr9/sYfFa1t/f5s1DdjeupXyajQQ5WMCj2MI6NCD1UyKvSgo9jiNI9R/jaPr42Is9vbjwGeB9wGfAB4Sfu0K4Dfb2+/u71P+/j72x3XGO3ZA73exrVer1mfhfpmTE/9KhkVeqiSUaGHcWRU6KFKRoUetAODJuxhbMA3ADcDtwC3Aj/drl8KfBC4A/htYK5dP7O9f0f7+KXHy/DI9Gh4ZfOpk1Ghh3FkVOihSkaFHsaRUaGHKhkVelCDSVyAOGpegDg6S0tLAOzdu3cm65sxPfWrZFTooUpGhR7GkVGhhyoZFXrQhC5AlCRJkipzmJYkSZI6cpiWJEmSOnKYliRJkjpymJYkSZI6cpiWJEmSOnKYliRJkjpymJYkSZI6cpiWJEmSOnKYliRJkjpymJYkSZI6cpiWJEmSOnKYliRJkjpymJYkSZI6cpiWJEmSOnKYliRJkrrKzJndFhYWUsO3vJw5N3c44ZHs95v7s1TfjOmpXyWjQg9VMir0MI6MCj1UyajQgxrASg6YRyc+EJ/M5jA9fMvLmb1e852xtvV6w/tgjrq+GdNTv0pGhR6qZFToYRwZFXqoklGhBx211TAdzWOzaXFxMVdWVia9G6XMz8PBg8euz83dx65du0+6/r5917C6esHI6psxPfWrZFTooUpGhR7GkVGhhyoZk+yh34cDB4YSoVZE7M/Mxc3rnjOtDQ4dGry+unr+UOpvVWdY9c2YnvpVMir0UCWjQg/jyKjQQ5WMSfaw1d/nGoFBh6tnZfM0j+Hr9zf+qGht6/dno74Z01O/SkaFHqpkVOhhHBkVeqiSUaEHHcUWp3l4ZFob7NkDvd7GtV6vWZ+F+mZMT/0qGRV6qJJRoYdxZFTooUpGhR60A4Mm7FnZPDI9Gl7ZfOpkVOhhHBkVeqiSUaGHcWRU6KFKRoUe1MALEHUilpaWANi7d+9M1jdjeupXyajQQ5WMCj2MI6NCD1UyKvQgL0CUJEmShs5hWpIkSerIYVqSJEnqyGFakiRJ6shhWpIkSerIYVqSJEnqyGFakiRJ6shhWpIkSerIYVqSJEnqyGFakiRJ6shhWpIkSerIYVqSJEnqyGFakiRJ6shhWpIkSerIYVqSJEnqyGFakiRJ6mjkw3REnBYRN0fE9e39p0TETRFxe0RcGxFntOtz7f072sfnR71vkiRJ0knJzJFuwL8G3gFc396/Dtjd3n4L8OPt7VcAb2lv7wauPV7thYWF1PAtL2fOzR1OeCT7/eb+LNU3Y3rqV8mo0EOVjAo9jCOjQg9VMir0oAawkoNm3UGLw9qAi4EbgW8DrgcCuB84vX382cD72tvvA57d3j69fV5sV99heviWlzN7veY7Y23r9Yb3wRx1fTOmp36VjAo9VMmo0MM4Mir0UCWjQg86aqthOprHRiMi3gX8LHAW8GrgB4F9mfnU9vFLgPdm5jMi4lbg+Zl5T/vYncCzMvP+reovLi7mysrKyPb/VDQ/DwcPHrs+N3cfu3btPun6+/Zdw+rqBSOrb8b01K+SUaGHKhkVehhHRoUeqmRMsod+Hw4cGEqEWhGxPzMXN6+P7JzpiHgh8JnM3L9+ecBTcwePra97ZUSsRMTKkSNHhrCnWu/QocHrq6vnD6X+VnWGVd+M6alfJaNCD1UyKvQwjowKPVTJmGQPW/19rhEYdLh6GBvNEel7gAPAfcAXgKvxNI+p1u9v/FHR2tbvz0Z9M6anfpWMCj1UyajQwzgyKvRQJaNCDzqKSZwz/dUQWOLoBYi/zcYLEF/R3n4lGy9AvO54dR2mh8/zx06djAo9jCOjQg9VMir0MI6MCj1UyajQg46apmH6UuCDwB3tYD3Xrp/Z3r+jffzS49V1mB4Nr2w+dTIq9DCOjAo9VMmo0MM4Mir0UCWjQg9qbDVMj/QCxFHzAsTRWVpaAmDv3r0zWd+M6alfJaNCD1UyKvQwjowKPVTJqNCDJnABoiRJklSdw7QkSZLUkcO0JEmS1JHDtCRJktSRw7QkSZLUkcO0JEmS1JHDtCRJktSRw7QkSZLUkcO0JEmS1JHDtCRJkuqvGUoAAAtqSURBVNSRw7QkSZLUkcO0JEmS1JHDtCRJktSRw7QkSZLUkcO0JEmS1JHDtCRJktSRw7QkSZLUVWbO7LawsJAavuXlzLm5wwmPZL/f3J+l+mZMT/0qGRV6qJJRoYdxZFTooUpGhR7UAFZywDw68YH4ZDaH6eFbXs7s9ZrvjLWt1xveB3PU9c2YnvpVMir0UCWjQg/jyKjQQ5WMCj3oqK2G6Wgem02Li4u5srIy6d0oZX4eDh48dn1u7j527dp90vX37buG1dULRlbfjOmpXyWjQg9VMir0MI6MCj1UyZhkD/0+HDgwlAi1ImJ/Zi5uXvecaW1w6NDg9dXV84dSf6s6w6pvxvTUr5JRoYcqGRV6GEdGhR6qZEyyh63+PtcIDDpcPSubp3kMX7+/8UdFa1u/Pxv1zZie+lUyKvRQJaNCD+PIqNBDlYwKPegotjjNwyPT2mDPHuj1Nq71es36LNQ3Y3rqV8mo0EOVjAo9jCOjQg9VMir0oB0YNGHPyuaR6dHwyuZTJ6NCD+PIqNBDlYwKPYwjo0IPVTIq9KAGXoCoE7G0tATA3r17Z7K+GdNTv0pGhR6qZFToYRwZFXqoklGhB3kBoiRJkjR0DtOSJElSRw7TkiRJUkcO05IkSVJHDtOSJElSRw7TkiRJUkcO05IkSVJHDtOSJElSRw7TkiRJUkcO05IkSVJHDtOSJElSRw7TkiRJUkcO05IkSVJHDtOSJElSRw7TkiRJUkcO05IkSVJXmTmyDTgAfAz4CLDSrp0L3ADc3v55TrsewJuAO4BbgMuOV39hYSE1fMvLmXNzhxMeyX6/uT9L9c2YnvpVMir0UCWjQg/jyKjQQ5WMCj2osTbLbt7GMUyft2ntjcBV7e2rgDe0ty8H3tsO1buAm45X32F6+JaXM3u95jtjbev1hvfBHHV9M6anfpWMCj1UyajQwzgyKvRQJaNCDzpqq2E6msdGIyIOAIuZef+6tU8CS5l5OCKeBOzNzK+PiF9rb79z8/O2qr+4uJgrKysj2/9T0fw8HDx47Prc3H3s2rX7pOvv23cNq6sXjKy+GdNTv0pGhR6qZFToYRwZFXqokjHJHvp9OHBgKBFqRcT+zFzcvD7qc6YT+KOI2B8RV7ZrT1wbkNs/z2/XLwLuXvfae9q1DSLiyohYiYiVI0eOjHDXT02HDg1eX109f/ADJ2irOsOqb8b01K+SUaGHKhkVehhHRoUeqmRMsoet/j7XCAw6XD2sDbiw/fN84KPAtwIPbXrOg+2ffwA8Z936jcDCdvU9zWP4+v2NPypa2/r92ahvxvTUr5JRoYcqGRV6GEdGhR6qZFToQUexxWkeIz0ynZn3tn9+Bvg94JnAp9vTO2j//Ez79HuAS9a9/GLg3lHun461Zw/0ehvXer1mfRbqmzE99atkVOihSkaFHsaRUaGHKhkVetAODJqwh7EBjwXOWnf7z4HnA/+JjRcgvrG9/Z1svADxg8fL8Mj0aCwvN/9HG5Eju7J5lPXNmJ76VTIq9FAlo0IP48io0EOVjAo9qMG4L0CMiEtpjkYDnA68IzP3RMQTgOuAJwOHgO/JzAciIoBfaQfuLwAvz8xtry70AkRJkiSNw1YXIJ4+qsDMvAv4xgHrfwk8d8B6Aq8c1f5IkiRJw+a/gChJkiR15DAtSZIkdeQwLUmSJHXkMC1JkiR15DAtSZIkdeQwLUmSJHXkMC1JkiR1NLJ/tGUcIuIIcHDS+1HYecD9k94JAb4X08L3YXr4XkwP34vp4XsxWv3M/NrNizM9TGu0ImJl0L/0o/HzvZgOvg/Tw/dievheTA/fi8nwNA9JkiSpI4dpSZIkqSOHaW3n1ye9A/oq34vp4PswPXwvpofvxfTwvZgAz5mWJEmSOvLItCRJktSRw7SOERHPj4hPRsQdEXHVpPenuoi4JCI+EBG3RcTHI+In2vVzI+KGiLi9/fOcdj0i4k3t+3NLRFw22Q7qiYjTIuLmiLi+vf+UiLipfS+ujYgz2vW59v4d7ePzk9zvaiLi7Ih4V0R8ov18PNvPxfhFxP/e/rfp1oh4Z0Sc6WdifCLi7RHxmYi4dd3aCX8OIuKK9vm3R8QVk+ilKodpbRARpwFvBl4APB34/oh4+mT3qrwvA/8mM58G7AJe2X7NrwJuzMyvA25s70Pz3nxdu10J/Or4d7m8nwBuW3f/DcAvtu/Fg8APt+s/DDyYmU8FfrF9nobnl4H/kZl/B/hGmvfEz8UYRcRFwKuAxcx8BnAasBs/E+P0m8DzN62d0OcgIs4FXgs8C3gm8Nq1AVwnz2Famz0TuCMz78rMvwGuAV404X0qLTMPZ+aH29ufoxkYLqL5uv9W+7TfAl7c3n4R8F+zsQ84OyKeNObdLisiLga+E3hrez+AbwPe1T5l83ux9h69C3hu+3ydpIj4GuBbgbcBZObfZOZD+LmYhNOBx0TE6UAPOIyfibHJzD8FHti0fKKfg38M3JCZD2Tmg8ANHDugqyOHaW12EXD3uvv3tGsag/ZHot8M3AQ8MTMPQzNwA+e3T/M9Gq1fAv4d8JX2/hOAhzLzy+399V/vr74X7eMPt8/XybsUOAL8l/aUm7dGxGPxczFWmfkXwM8Dh2iG6IeB/fiZmLQT/Rz4+Rghh2ltNugIgr/yZQwi4nHA7wD/KjM/u91TB6z5Hg1BRLwQ+Exm7l+/POCpuYPHdHJOBy4DfjUzvxn4PEd/lD2I78UItKcCvAh4CnAh8FiaUwk28zMxHbb6+vu+jJDDtDa7B7hk3f2LgXsntC+njIh4NM0gfXVm/m67/Om1H1O3f36mXfc9Gp1/AHx3RBygOcXp22iOVJ/d/ogbNn69v/petI8/nmN/HKtu7gHuycyb2vvvohmu/VyM1/OAT2Xmkcz8EvC7wN/Hz8SknejnwM/HCDlMa7MPAV/XXql9Bs2FJu+e8D6V1p5P+Dbgtsz8hXUPvRtYu+L6CuD3162/rL1qexfw8NqP+3RyMvMnM/PizJyn+d5/f2b+APAB4CXt0za/F2vv0Uva53u0Zwgy8z7g7oj4+nbpucD/i5+LcTsE7IqIXvvfqrX3wc/EZJ3o5+B9wHdExDntTxu+o13TEPiPtugYEXE5zdG404C3Z+aeCe9SaRHxHOB/Ah/j6Hm6/57mvOnrgCfT/IX2PZn5QPsX2q/QXDzyBeDlmbky9h0vLiKWgFdn5gsj4lKaI9XnAjcD/ywzVyPiTOC/0Zzn/gCwOzPvmtQ+VxMR30RzIegZwF3Ay2kOAvm5GKOIeD3wfTS/eehm4Edozrf1MzEGEfFOYAk4D/g0zW/l+O+c4OcgIn6I5u8WgD2Z+V/G2UdlDtOSJElSR57mIUmSJHXkMC1JkiR15DAtSZIkdeQwLUmSJHXkMC1JkiR15DAtSZIkdeQwLUnFRcR8RHwxIj5ygq/7voi4IyKuH9W+SdKsc5iWpFPDnZn5TSfygsy8luYf6JAkbcFhWpJmWER8S0TcEhFnRsRjI+LjEfGM47xmPiI+ERFvjYhbI+LqiHheRPyviLg9Ip45rv2XpFl3+qR3QJLUXWZ+KCLeDfyfwGOA5cy8dQcvfSrwPcCVwIeAlwLPAb6b5p8cfvFo9liSanGYlqTZ9zM0A/FfA6/a4Ws+lZkfA4iIjwM3ZmZGxMeA+ZHspSQV5GkekjT7zgUeB5wFnLnD16yuu/2Vdfe/ggdaJGnHHKYlafb9OvBTwNXAGya8L5J0SvHogyTNsIh4GfDlzHxHRJwG/HlEfFtmvn/S+yZJp4LIzEnvgyRphCJiHrg+M7f9LR9bvHYJeHVmvnDIuyVJJXiahyTV9wjw+C7/aAvwn4EHR7JXklSAR6YlSZKkjjwyLUmSJHXkMC1JkiR15DAtSZIkdeQwLUmSJHXkMC1JkiR19P8DpbkyF4oMehsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Define figure size\n",
    "rcParams['figure.figsize'] = 12, 5\n",
    "\n",
    "# Calculate coordinates of grid points\n",
    "NX, NZ, x, z, XX, ZZ = coord_def(Xmax,Zmax,dh)\n",
    "\n",
    "# Plot Cartesian mesh\n",
    "plt.plot(XX, ZZ, 'k')\n",
    "plt.plot(XX.T, ZZ.T, 'k')\n",
    "\n",
    "# Plot grid points\n",
    "plt.plot(XX, ZZ, 'bo')\n",
    "\n",
    "plt.title(\"Equidistant 2D Cartesian grid\" )\n",
    "plt.xlabel(\"x [m]\")\n",
    "plt.ylabel(\"z [m]\")\n",
    "plt.axis('equal')\n",
    "plt.gca().invert_yaxis()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nice, next we create and visualize a homogeneous Vp-model discretized on this Cartesian grid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArQAAAFNCAYAAADxS/yWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZhlVX3v//eHZlJAGRoVaCZvSJSYBAwiibmKaBSIARQnMgCGXDJoEqMYwMcrSuSq9zEhGo2GKwgoMgTlZ8dgCJdBrlEQUGRUaUGlbRQJk4qg0N/fH3sVHIrqrnOq63Tt6nq/fPZTe689nHX2qYOfXrX2WqkqJEmSpPlqvbmugCRJkrQmDLSSJEma1wy0kiRJmtcMtJIkSZrXDLSSJEma1wy0kiRJmtcMtJI0jyS5NMkfr+E1PpLkf85WncYlyd5Jlg957DuSfGLcdZLUTwZaaZYl+XaSF08qOzzJF+aqTuuiFmB+nuTHSe5J8sUkvzHX9ZoPqupPq+pvYbTQKEl9ZaCVNJ+dXVWbAlsDXwA+nSRzXCdJ0lpmoJXmQJJntj8d35PkhiQHDOw7Nck/Jflca338zyRPS/IPSe5O8vUkuw95ra2S/GuS+5JcmeRdgy3FSZ6R5MIkdyX5RpJXT6rHh5L8W5IfJbkiyX8b8twnJzk9yQ+TfCfJ25Ks1/Y95k/DSXZKUknWb9uHJ7mlveatSX5/uvtZVT8HTgOeBmw1xf1+Z5J/bOsbJPlJkv/dtp+Q5IEkW7Ttf0ny/ST3JrksyS+38r1a+aKB6748ybVtfb0kxyT5VpL/SnJOki2nqm+Sm5K8bGB7/SR3Jnn2wGt9sX2mX0uy9yqus167t99Jcke7508e2P9bA9e5LcnhrfzU9ruwCfA5YNv2u/bjJNsmuT/JVgPX+fX2WW4wRR3e0e7ZJ9pndl2SX0xybKvTbUleMnD8tkmWtt+bZUn+x8C+J7S63Z3kRuA5k15r2ySfanW5NclfTnVfJC08BlppLWuh4F+B/wCeAvwFcEaSXxo47NXA24DFwIPAl4CvtO1zgb8f8lofAn5CF/QOa8tEPTYBLgQ+2c49BPiniQDXHAK8E9gCWAacMOS5/wg8GXg68ALgUOB1Q9ybTYAPAPtV1WbAbwLXDHHeRsDhwPKqunOKQz4P7N3WnwN8v9UL4DeAb1TV3W37c8Au7X19BTgDoKoup7uX+wxc9/fo7gHAXwIHtetuC9xNd/+ncibdPZvwUuDOqvpKku2AfwPeBWwJHAV8KsnWU1zn8La8kO5ebwp8ECDJDu29/CNdC/ZuTLqXVfUTYD9gRVVt2pYVwKV0v4MT/gA4q/3DYSq/C3yc7vfkq8AFdP//sh1wPPDPk977crp79ErgfyV5Udt3HPDf2vJSHvv7uh7d7/rX2nVfBLwxyUtXUSdJC0lVubi4zOICfBv4MXDPwHI/8IW2/7/TBar1Bs45E3hHWz8V+D8D+/4CuGlg+1eAe6a7FrAI+DnwSwP73jVQj9cA/29S3f8ZOG6gHh8d2Lc/8PXpzm2v+yCw68C+PwEubevvAD4xsG8noID1gU3a/ToYeMI09/kdwM/a8XcAFwO/vopjnwA8QNd6ewzwVrpQtSldYP/AKs7bvNXtyQP375S2vhldwN2xbd8EvGjg3G3a/V9/iuv+AvAj4Ilt+wzg7W39aODjk46/ADisrV8K/HFbvwj484HjfmniNYFjgfNW8b5OBd7V1vem+4fA4P7XAP/Z1he137E9V/M5XDiw/bt0v/+LBu5TtXu5PfAwsNnA8e8GTm3rtwD7Duw7cqJuwHOB70567WOBj031e+Xi4rKwFltopfE4qKo2n1iAPx/Yty1wW1WtHCj7Dl2r04QfDKz/dIrtTYe41tZ0wea2gX2D6zsCz21/jr4nyT3A79O15k74/sD6/QOvu7pzFwMbtnqs6v1NqboWw9cAfwrcnq67wzNWc8o57R4/par2qaqrAdJ1vZj4E/p/r6qfAlfRtZ4+n67F9ovA81rZ59t5i5K8p3UbuI/uHye09wRda+wrWovwK4CvVNXE+9wROG/gftxEF96eOsX7XNb2/26SJwIH8GhL747Aqybd29+iC8iTbcvj7/P67TW3B761mnu3Op8Bdk3ydOC3gXur6surOX7y7+edVfXwwDZ0vzvbAndV1Y8m1Xnid2NbHvs7OvjedqTrGjF4X97KFPdX0sKz/lxXQFqAVgDbJ1lvIIjuAHxzlq/1Q+AhYMnAtbcfOPc24PNV9dszeN1Vntv6mP6cLoDcOFCn77X1nwBPHDhlMEBTVRcAFyR5Al2L6P+ha4keWlX98hTFn6frLrA7cGXbfimwJ3BZO+b3gAOBF9OF2SfTdR1Iu+6NSb5D92f6we4G0N2TP6qq/xyymhPdDtYDbmwhd+I6H6+q/7HKMx+1gu4+T9iB7jP/QbvOnkNcox5XUPVAknPo/pHyDLruBLNhBbBlks0GQu3g78btdL+jNwzsm3AbcGtV7TJLdZG0DrGFVlr7rqALdX+T7gGlven+THvWbF6rtZB9GnhHkie2ls5DB879LPCLSf6wnbtBkuckeeYQr7vKc9vrngOckGSzJDsCbwImHgS7Bnh+kh3aA0zHTlw0yVOTHND60j5I96frh5kdn6d7/zdW1c9of7qnC0k/bMds1l73v+hC9/+a4jqfpOsv+3zgXwbKP0L3nnds72XrJAeupj5nAS8B/ozHBuNP0LXcvrS1GG+cbmitJVNc40zgr5PsnGTTVt+zq+ohum4ML07y6nQPnW2VZLcprvEDYKsMPEzWnE7XP/cAHv3s1khV3UbXMv7u9r5+FTii1RW635tjk2zR3u9fDJz+ZeC+JEe3h8cWJXlWksc8OCZpYTLQSmtZC1MH0LXy3Qn8E3BoVX19DNd6A10r4/fpWtnOpAtstBaylwCvpWs5+z7wXmCjIV53unP/gi5o30I3nNYngVPauRcCZwPXAlfTheMJ6wFvbte8i647wGB3jTXxRbq+tBOtsTfS9au9bOCY0+n+zP29tv/yKa5zJl2/04vrsQ+gvR9YCvxHkh+1c5+7qspU1e10D/v9Jt39mCi/ja6V+K10rey3AW9h6v9en0L3uV4G3Nrez1+063yXrt/zm+nu5TXAr01Rj6+393RL+1P+tq38P4GVdN0qvr2q9zEDh9D1m14BnEfXZ/vCtu+ddPf/VroHHR9pGW7/UPpduofbbqX7ff8o3e+3pAUuVY/7a5OkdVSS9wJPq6rDpj1YC16Si4FPVtVH57oukrQ6ttBK67B0Y8X+ajp70v1597y5rpf6r/0p/9kMtB5LUl/1KtAm2TfdAO3Lkhwz1/WR1gGb0fWj/Qld/8S/o3uCXVqlJKcB/xd446QRCSTNc0lOSTfpyfWr2P+MJF9K8mCSoybtmzKntX78VyS5OcnZSTYc9/t4XL370uWgPRn9TbohYpbTPYV8SFXduNoTJUmSNJQkz6d74Pb0qnrWFPufQjd6ykHA3VX1vla+ypzWRkX5dFWdleQjwNeq6sNr5x11+tRCuyewrKpuaQ+6nEX3YIQkSZJmQVVdRveg6Kr231FVV9INvzhoypyWJHRDIp7bjjuNLgyvVX0KtNvx2AG1lzPEQOySJEkau1XltK3oZq98aFL5WtWniRUyRdnj+kMkOZJuOkSy4Ya/vsFTnzLuekmSJM3Yz25bfmdVbT3KOS994Sb1X3eNPgz31dc+eAPdEH4TTqqqk0a+0OOtKqcNld/GrU+BdjmPncVoCd04hY/RPpSTADbaYfva7s1/vXZqJ0mSNAO3vvHN35n+qMf6r7se5ssX7DD9gZMs2ubmB6pqj5FPnN6qctqdwOZJ1m+ttFPmt3HrU5eDK4Fd2pNyG9IN2L50juskSZK01hWwcgb/G6Mpc1p1owtcAryyHXcYczCaTm9aaKvqoSRvAC4AFgGnVNUN05wmSZK0DioertkPqEkmZjtcnGQ5cBywAUBVfSTJ04CrgCcBK5O8Edi1qu5bTU47GjgrybuArwInz3rFp9GbQAtQVecD5891PSRJkuZS10I7+11Rq+qQafZ/n67bwFT7psxpVXUL3SgIc6ZXgVaSJEmdMXchWKcYaCVJknqmKB7uyeRX84GBVpIkqYfG0eVgXWWglSRJ6pkCHjbQDs1AK0mS1EO20A7PQCtJktQzBfahHYGBVpIkqYcc42B4BlpJkqSeKco+tCMw0EqSJPVNwcPm2aEZaCVJknqmmylMwzLQSpIk9U54mMx1JeYNA60kSVLPFLDSLgdDM9BKkiT1kC20w1tvrisgSZIkrQlbaCVJknqmm/rWFtphGWglSZJ6aGUZaIdloJUkSeoZW2hHY6CVJEnqmSI87KNOQzPQSpIk9ZBdDoZnoJUkSeoZuxyMxkArSZLUO+HhssvBsAy0kiRJPVPASvvQDs1AK0mS1EN2ORiegVaSJKlnquxyMAoDrSRJUg+ttIV2aEZ/SZKknulGOVhv5GU6SU5JckeS61exP0k+kGRZkmuTPLuVvzDJNQPLA0kOavtOTXLrwL7dZvNeDMMWWkmSpN4ZW5eDU4EPAqevYv9+wC5teS7wYeC5VXUJsBtAki2BZcB/DJz3lqo6dxwVHoaBVpIkqWfGNcpBVV2WZKfVHHIgcHpVFXB5ks2TbFNVtw8c80rgc1V1/6xXcIbsciBJktRDD1dGXmbBdsBtA9vLW9mg1wJnTio7oXVRODHJRrNRkVEYaCVJknqmyEz70C5OctXAcuSILz1VKq5HdibbAL8CXDCw/1jgGcBzgC2Bo0d8zTVmlwNJkqR1x51VtccanL8c2H5gewmwYmD71cB5VfXziYKB7ggPJvkYcNQavP6M2EIrSZLUQytrvZGXWbAUOLSNdrAXcO+k/rOHMKm7QWu1JUmAg4ApR1AYJ1toJUmSemZi2K7ZluRMYG+6rgnLgeOADQCq6iPA+cD+dKMY3A+8buDcnehabz8/6bJnJNmarrvCNcCfznrFp2GglSRJ6pli1h7yeux1qw6ZZn8Br1/Fvm/z+AfEqKp9ZqVya8BAK0mS1EPjGLZrXWWglSRJ6pkqxjWxwjrJQCtJktQ7YeWUI2hpKgZaSZKknilsoR2FgVaSJKmHxjHKwbrKQCtJktQzRVg5hlEO1lUGWkmSpB6yhXZ4BlpJkqSeKZitmb8WBAOtJElS74SHHeVgaAZaSZKknrGFdjQGWkmSpB6yhXZ4BlpJkqSeqYottCMY251KckqSO5JcP1C2ZZILk9zcfm7RypPkA0mWJbk2ybPHVS9JkqT54OFab+RloRrnOz8V2HdS2THARVW1C3BR2wbYD9ilLUcCHx5jvSRJkrQOGVugrarLgLsmFR8InNbWTwMOGig/vTqXA5sn2WZcdZMkSeqzAlaSkZeFam33oX1qVd0OUFW3J3lKK98OuG3guOWt7PbJF0hyJF0rLou22GK8tZUkSZoTWdBdCEbVl4fCpvonRU11YFWdBJwEsNEO2095jCRJ0nzWDdu1cFtcR7W2A+0PkmzTWme3Ae5o5cuB7QeOWwKsWMt1kyRJ6g2nvh3e2r5TS4HD2vphwGcGyg9tox3sBdw70TVBkiRpoSnCyhp9WajG1kKb5Exgb2BxkuXAccB7gHOSHAF8F3hVO/x8YH9gGXA/8Lpx1UuSJGk+WGkL7dDGFmir6pBV7HrRFMcW8Ppx1UWSJGk+qYKHF3CL66j68lCYJEmSBizkLgSjMtBKkiT1TNeH1i4HwzLQSpIk9dDDC3iihFEZaCVJknrGcWhHY1u2JElS73RdDkZdpr1qckqSO5Jcv4r9SfKBJMuSXJvk2QP7Hk5yTVuWDpTvnOSKJDcnOTvJhrNyC0ZgoJUkSeqhlWTkZQinAvuuZv9+wC5tORL48MC+n1bVbm05YKD8vcCJVbULcDdwxCjvczYYaCVJknpmYtiuUZfpr1uXAXet5pADgdOrczmweZvddUpJAuwDnNuKTgMOGvqNzhIDrSRJUg+No8vBELYDbhvYXt7KADZOclWSy5NMhNatgHuq6qEpjl9rfChMkiSpZyamvp2BxUmuGtg+qapOGuH8qV602s8dqmpFkqcDFye5DrhvNcevNQZaSZKkdcedVbXHGpy/HNh+YHsJsAKgqiZ+3pLkUmB34FN03RLWb620jxy/NtnlQJIkqYfG9FDYdJYCh7bRDvYC7q2q25NskWQjgCSLgecBN1ZVAZcAr2znHwZ8ZjYqMgpbaCVJknpmXOPQJjkT2Juua8Jy4DhgA4Cq+ghwPrA/sAy4H3hdO/WZwD8nWUnXIPqeqrqx7TsaOCvJu4CvAifPesWnYaCVJEnqoXFMfVtVh0yzv4DXT1H+ReBXVnHOLcCes1LBGTLQSpIk9U3N+KGwBclAK0mS1DMFs9UndkEw0EqSJPWQLbTDM9BKkiT1zLgeCltXGWglSZJ6yEA7PAOtJElSz6zBTGELkoFWkiSph3wobHgGWkmSpL4puxyMwkArSZLUMz4UNhoDrSRJUg8ZaIdnoJUkSeoZHwobjYFWkiSph8pAOzQDrSRJUg85ysHwDLSSJEk9U45yMJL15roCkiRJ0pqwhVaSJKmH7EM7PAOtJElS7zjKwSgMtJIkST1kC+3wDLSSJEk940xhozHQSpIk9U11Ix1oOAZaSZKkHnIc2uEZaCVJknqmsA/tKAy0kiRJveMoB6Mw0EqSJPWQfWiHZ6CVJEnqIbscDM+pbyVJknqmqgu0oy7TSXJKkjuSXL+K/UnygSTLklyb5NmtfLckX0pyQyt/zcA5pya5Nck1bdlt1m7EkGyhlSRJ6qEx9aE9FfggcPoq9u8H7NKW5wIfbj/vBw6tqpuTbAtcneSCqrqnnfeWqjp3HBUehoFWkiSph8bRh7aqLkuy02oOORA4vaoKuDzJ5km2qapvDlxjRZI7gK2Be1Z1obXJLgeSJEk9NI4uB0PYDrhtYHt5K3tEkj2BDYFvDRSf0LoinJhko9moyCgMtJIkST1TjB5mW6BdnOSqgeXIEV96qlT8SFtxkm2AjwOvq6qVrfhY4BnAc4AtgaNHfsNryC4HkiRJPTTDHgd3VtUea/Cyy4HtB7aXACsAkjwJ+DfgbVV1+cQBVXV7W30wyceAo9bg9WfEFlpJkiRNWAoc2kY72Au4t6puT7IhcB5d/9p/GTyhtdqSJMBBwJQjKIzT2Fpok2xP9wTd04CVwElV9f4kWwJnAzsB3wZeXVV3t5vwfmB/uifpDq+qr4yrfpIkSb1V4xmHNsmZwN50XROWA8cBGwBU1UeA8+my2DK6PPa6duqrgecDWyU5vJUdXlXXAGck2Zquu8I1wJ/OesWnMc4uBw8Bb66qryTZjG54hwuBw4GLquo9SY4BjqHra7GqYSIkSZIWnvGMcnDINPsLeP0U5Z8APrGKc/aZndrN3Ni6HFTV7RMtrFX1I+AmuqfkDgROa4edRtc0DQPDRLR+GZtPNGFLkiQtNHM0ysG8tFb60LbxznYHrgCeOtF5uP18Sjts2mEiJEmSFoputrDRloVq7KMcJNkU+BTwxqq6r+sqO/WhU5Q97qNpw08cCbBoiy1mq5qSJEm9UYynD+26aqwttEk2oAuzZ1TVp1vxDwaehtsGuKOVr3KYiEFVdVJV7VFVeyzadJPxVV6SJGmuFFAZfVmgxhZo26gFJwM3VdXfD+xaChzW1g8DPjNQ/rhhIsZVP0mSpD6zy8Hwxtnl4HnAHwLXJbmmlb0VeA9wTpIjgO8Cr2r7VjVMhCRJ0sKzgAPqqMYWaKvqC0zdLxbgRVMcP+UwEZIkSQvPwh61YFROfStJktRHttAOzUArSZLUN2OaKWxdZaCVJEnqowXQQpvkA0Mcdl9VvW11BxhoJUmSemlBtNAeCLx9mmOOAQy0kiRJ884CaKEFTqyq01Z3QJJpZ9Iy0EqSJPXRAgi0VfUPs3HMWGcKkyRJ0gwssJnCkvzvJE9KskGSi5LcmeQPhj3fQCtJkqS59pKqug94GbAc+EXgLcOebJcDSZKkHlpgU9lu0H7uD5xZVXclw7c4G2glSZL6aGEF2n9N8nXgp8CfJ9kaeGDYk+1yIEmS1EcLoA9tkm0AquoY4DeAParq58D9dEN6DWW1LbRJ7puuHsDtVfWLw76gJEmSppeF0UJ7ShuW61Lg34EvAFTVT4CfDHuR6bocfKuqdl/dAUm+OuyLSZIkaQjFguhyUFX7JdkY2Bt4OfC+JN+lC7f/XlXfHeY60wXag4e4xjDHSJIkaWjzswvBTFTVA7QAC5BkZ2A/4INJnlZVe053jdUG2qq6ZXA7yZMGz6mquyYfI0mSpFmwAFpoJ2tZ817grLb8eJjzhhrlIMmfAMfTPXk2cXsLePrINZUkSdL0FlCgXVXWrKqhsuaww3YdBfxyVd05ehUlSZI0sgUUaFnDrDlsoP0W3fAJkiRJGreJqW8XjjXKmsMG2mOBLya5AnhworCq/nKmLyxJkqRVG8ewXUlOoZte9o6qetYU+wO8n27GrvuBw6vqK23fYcDb2qHvqqrTWvmvA6cCTwDOB/6qauR5ztYoaw4baP8ZuBi4Dlg5YgUlSZI0qvF0OTgV+CBw+ir27wfs0pbnAh8GnptkS+A4YI9Ws6uTLK2qu9sxRwKX0wXafYHPjVivNcqawwbah6rqTaNeXJIkSf1RVZcl2Wk1hxwInN5aWC9PsnmbzWtv4MKqugsgyYXAvkkuBZ5UVV9q5acDBzF6oF2jrDns1LeXJDkyyTZJtpxYZvqikiRJWr3U6Mss2A64bWB7eStbXfnyKcpHtUZZc9gW2t9rP48dKHPYLkmSpHGZ2UNhi5NcNbB9UlWdNML5U71ozaB8VGuUNYcKtFW184iVkiRJ0kzNfOrbO6tqjzV45eXA9gPbS4AVrXzvSeWXtvIlUxw/kjXNmqvtcpDk2dNdYJhjJEmSNC8sBQ5NZy/g3qq6HbgAeEmSLZJsAbwEuKDt+1GSvdoICYcCnxn2xWYra07XQvuxJHszdXPyhJOB3ad7IUmSJI1gPMN2nUnX0ro4yXK6kQs2AKiqj9CNUrA/sIxu2K7XtX13Jflb4Mp2qeMnHhAD/oxHh+36HKM9EDYrWXO6QPtk4OppXuSH01xDkiRJIxrHOLRVdcg0+wt4/Sr2nQKcMkX5VcDjxrQd0qxkzdUG2qraabQ6SZIkaVYsgKlvZytrDjvKgSRJktamBRBoZ4uBVpIkqWdmcVzZBcFAK0mS1EczG4d2QRpqprAkFyXZf1LZKIP0SpIkaRQ1g2UeS/KKJH+f5O+SvHyUc4ed+nZn4Ogkxw2UrcmgvZIkSVqNOZr6dk4k+SfgT4HrgOuBP0nyoWHPH7bLwT3Ai4APJPlX4A9GragkSZJGMI8D6gy8AHhWGzaMJKfRhduhDNtCm6p6qKr+HPgU8AXgKaPWVJIkSUOYQevsfG6hBb4B7DCwvT1w7bAnD9tC+5GJlao6Ncl1rGLQXUmSJM2C+R1QR7UVcFOSL7ft5wBfSrIUoKoOWN3JQwXaqvrnSdtXA380el0lSZI0lIUVaN++Jic7bJckSVIPzfMuBENJ8kHgk1X1+TW5zrB9aCVJkqTZdjPwd0m+neS9SXabyUUMtJIkSX20AMahrar3V9Vv0I1ycBfwsSQ3JXl7kl8c9joGWkmSpL5ZYKMcVNV3quq9VbU78HvAy4Gbhj3fQCtJkqQ5lWSDJL+b5Azgc8A3gYOHPd+HwiRJkvpoHre4DivJbwOHAL8DfBk4Cziyqn4yynUMtJIkSX20AAIt8Fbgk8BRVXXXTC9ioJUkSeqZML/7xA6rql44G9cZWx/aJBsn+XKSryW5Ick7W/nOSa5IcnOSs5Ns2Mo3atvL2v6dxlU3SZKk3lsAoxzMlnE+FPYgsE9V/RqwG7Bvkr2A9wInVtUuwN3AEe34I4C7q+oXgBPbcZIkSQvPAhvlYE2NLdBW58dtc4O2FLAPcG4rPw04qK0f2LZp+1+UJOOqnyRJUq/ZQju0sQ7blWRRkmuAO4ALgW8B91TVQ+2Q5cB2bX074DaAtv9eYKtx1k+SJKm3DLRDG+tDYVX1MLBbks2B84BnTnVY+zlVa+zjPpokRwJHAizaYotZqqkkSVK/LOQuBKNaKxMrVNU9wKXAXsDmSSaC9BJgRVtfDmwP0PY/mW4KtMnXOqmq9qiqPRZtusm4qy5JkjQ3bKEd2jhHOdi6tcyS5AnAi+mmMLsEeGU77DDgM219adum7b+4qhbwRyNJkhasmYTZBZyaxtlCuw1wSZJrgSuBC6vqs8DRwJuSLKPrI3tyO/5kYKtW/ibgmDHWTZIkqdfGMcpBkn2TfKMNk/q4rJVkxyQXJbk2yaVJlrTyFya5ZmB5IMlBbd+pSW4d2LfbbN+L6YytD21VXQvsPkX5LcCeU5Q/ALxqXPWRJEmaV2a5xTXJIuBDwG/TdfW8MsnSqrpx4LD3AadX1WlJ9gHeDfxhVV1CNwwrSbYElgH/MXDeW6rqXObIWulDK0mSpNGMoYV2T2BZVd1SVT8DzqIbNnXQrsBFbf2SKfZD1zX0c1V1/8zf3ewy0EqSJPXR7PehfWSI1GZw+NQJXwMObusvBzZLMnkY1dcCZ04qO6F1UzgxyUbT1mSWGWglSZL6ZuYPhS1OctXAcuTAVYcZIvUo4AVJvgq8APgeMDF/AEm2AX4FuGDgnGOBZwDPAbake15qrRrrOLSSJEkaXZg6fQ7hzqraYxX7HhkitRkcPhWAqloBvAIgyabAwVV178AhrwbOq6qfD5xze1t9MMnH6ELxWmULrSRJUh/NfpeDK4FdkuycZEO6rgNLBw9IsjjJRD48Fjhl0jUOYVJ3g9ZqS5IABwHXD/sWZ4uBVpIkaQGoqoeAN9B1F7gJOKeqbkhyfJID2mF7A99I8k3gqcAJE+cn2Ymuhffzky59RpLrgOuAxcC7xvg2pmSXA0mSpB4ax9S3VXU+cP6ksrcPrJ8LTDn8VlV9m8c/REZV7TO7tRydgVaSJKmPFvDMX6My0EqSJPWRgXZoBlpJkqS+GXIqW3UMtJIkSX1koB2agVaSJKmHbKEdnoFWkiSpjwy0QzPQSpIk9ZAttMMz0EqSJPXNcDN/qTHQSlOiMDQAABAMSURBVJIk9ZGBdmgGWkmSpJ4JdjkYhYFWkiSpjwy0QzPQSpIk9VDKRDssA60kSVLf+FDYSAy0kiRJPWQf2uEZaCVJkvrIQDu09ea6ApIkSdKasIVWkiSph+xyMDwDrSRJUh8ZaIdmoJUkSeqbsoV2FAZaSZKkPjLQDs1AK0mS1DNOfTsaA60kSVIfOVPY0Ay0kiRJPWQL7fAMtJIkSX3j1LcjMdBKkiT1UFbOdQ3mD2cKkyRJ6qOawTKNJPsm+UaSZUmOmWL/jkkuSnJtkkuTLBnY93CSa9qydKB85yRXJLk5ydlJNlyTtz0TBlpJkqQeSo2+rPZ6ySLgQ8B+wK7AIUl2nXTY+4DTq+pXgeOBdw/s+2lV7daWAwbK3wucWFW7AHcDR6zRG58BA60kSVLfFN0oB6Muq7cnsKyqbqmqnwFnAQdOOmZX4KK2fskU+x8jSYB9gHNb0WnAQcO/0dlhoJUkSeqh2W6hBbYDbhvYXt7KBn0NOLitvxzYLMlWbXvjJFcluTzJRGjdCrinqh5azTXHzofCJEmS+mhmoxwsTnLVwPZJVXVSW88Qr3IU8MEkhwOXAd8DJsLqDlW1IsnTgYuTXAfcN2s1XwMGWkmSpJ5Zg5nC7qyqPVaxbzmw/cD2EmDF4AFVtQJ4BUCSTYGDq+regX1U1S1JLgV2Bz4FbJ5k/dZK+7hrrg12OZAkSeqbmfSfnb4P7ZXALm1Ugg2B1wJLBw9IsjjJRD48FjillW+RZKOJY4DnATdWVdH1tX1lO+cw4DOzcAdGYqCVJElaAFoL6huAC4CbgHOq6oYkxyeZGLVgb+AbSb4JPBU4oZU/E7gqydfoAux7qurGtu9o4E1JltH1qT15rbyhAXY5kCRJ6qFxTH1bVecD508qe/vA+rk8OmLB4DFfBH5lFde8hW4EhTljoJUkSeojp74dmoFWkiSph8bRQruuMtBKkiT1TQErTbTDMtBKkiT1kXl2aAZaSZKkHrLLwfDGPmxXkkVJvprks2175yRXJLk5ydltHDSSbNS2l7X9O427bpIkSb01++PQrrPWxji0f0U31tmE9wInVtUuwN3AEa38CODuqvoF4MR2nCRJ0oKUGn1ZqMYaaJMsAX4H+GjbDrAPj45vdhpwUFs/sG3T9r+oHS9JkrSw1AyXBWrcLbT/APwNsLJtbwXc02aqgG5O4e3a+nbAbfDITBb3tuMfI8mRSa5KctXDP/7JOOsuSZI0JwKkauRloRpboE3yMuCOqrp6sHiKQ2uIfY8WVJ1UVXtU1R6LNt1kFmoqSZLUQytnsCxQ4xzl4HnAAUn2BzYGnkTXYrt5kvVbK+wSYEU7fjmwPbA8yfrAk4G7xlg/SZKk3lrILa6jGlsLbVUdW1VLqmon4LXAxVX1+8AlwCvbYYcBn2nrS9s2bf/FVX6SkiRpAbIP7UjWxigHkx0NvCnJMro+sie38pOBrVr5m4Bj5qBukiRJPTCDIbsWcDvgWplYoaouBS5t67cAe05xzAPAq9ZGfSRJkvpuIQ/DNSpnCpMkSeqjBdziOqq56HIgSZIkzRpbaCVJkvqmIAt4GK5RGWglSZL6yC4HQzPQSpIk9ZF5dmgGWkmSpB5yYoXhGWglSZL6yEA7NAOtJElS3xTgQ2FDM9BKkiT1TCi7HIzAQCtJktRHBtqhGWglSZL6yEA7NGcKkyRJ6puJPrSjLtNIsm+SbyRZluSYKfbvmOSiJNcmuTTJkla+W5IvJbmh7XvNwDmnJrk1yTVt2W0N3/3IbKGVJEnqodnuQ5tkEfAh4LeB5cCVSZZW1Y0Dh70POL2qTkuyD/Bu4A+B+4FDq+rmJNsCVye5oKruaee9parOndUKj8AWWkmSpD6qGn1ZvT2BZVV1S1X9DDgLOHDSMbsCF7X1Syb2V9U3q+rmtr4CuAPYepbe6Roz0EqSJPXODMLs9IF2O+C2ge3lrWzQ14CD2/rLgc2SbDV4QJI9gQ2Bbw0Un9C6IpyYZKNR3+2aMtBKkiT1TTHTQLs4yVUDy5EDV80qXmnQUcALknwVeAHwPeChRy6QbAN8HHhdVU302j0WeAbwHGBL4OhZuAMjsQ+tJElSH81sYoU7q2qPVexbDmw/sL0EWDF4QOtO8AqAJJsCB1fVvW37ScC/AW+rqssHzrm9rT6Y5GN0oXitsoVWkiSph1I18jKNK4FdkuycZEPgtcDSx7xmsjjJRD48FjillW8InEf3wNi/TDpnm/YzwEHA9Wv41kdmoJUkSVoAquoh4A3ABcBNwDlVdUOS45Mc0A7bG/hGkm8CTwVOaOWvBp4PHD7F8FxnJLkOuA5YDLxr7byjR9nlQJIkqY/GMLFCVZ0PnD+p7O0D6+cCjxt+q6o+AXxiFdfcZ5arOTIDrSRJUt8UsNKZwoZloJUkSeqdoYbhUmOglSRJ6iMD7dAMtJIkSX1koB2agVaSJKlv7EM7EgOtJElS7xTUzGZWWIgMtJIkSX1kl4OhGWglSZL6xi4HIzHQSpIk9ZEttEMz0EqSJPWRgXZoBlpJkqTecWKFURhoJUmS+qaAlY5yMCwDrSRJUh/ZQjs0A60kSVIfGWiHZqCVJEnqnXLYrhEYaCVJkvqmoJwpbGjrzXUFJEmSpDVhC60kSVIf2eVgaAZaSZKkPvKhsKEZaCVJkvqmynFoR2CglSRJ6iNbaIdmoJUkSeqhsoV2aAZaSZKk3ilbaEdgoJUkSeqbwlEORmCglSRJ6iMnVhjaWCdWSPLtJNcluSbJVa1syyQXJrm5/dyilSfJB5IsS3JtkmePs26SJEl9VUCtrJGX6STZN8k3Wt46Zor9Oya5qGWxS5MsGdh3WMtvNyc5bKD811veW9ayXGbrPgxrbcwU9sKq2q2q9mjbxwAXVdUuwEVtG2A/YJe2HAl8eC3UTZIkqX+quhbaUZfVSLII+BBd5toVOCTJrpMOex9welX9KnA88O527pbAccBzgT2B4yYaJeky25E8muP2nY1bMIq5mPr2QOC0tn4acNBA+enVuRzYPMk2c1A/SZKkOTeGFto9gWVVdUtV/Qw4iy5/DdqVrsER4JKB/S8FLqyqu6rqbuBCYN+W1Z5UVV+qqgJO59Fst9aMO9AW8B9Jrk5yZCt7alXdDtB+PqWVbwfcNnDu8lYmSZK08MxyCy3DZa2vAQe39ZcDmyXZajXnbtfWV3fNsRv3Q2HPq6oVSZ4CXJjk66s5dqr+Fo/7p0YLxhPh+MFb3/jm62ehnlr7FgN3znUlNCN+dvObn9/85Wc3f/3SqCf8iLsv+L917uIZvNbGE88tNSdV1UltfZisdRTwwSSHA5cB3wMeWs25Q+W3cRtroK2qFe3nHUnOo2vq/kGSbarq9tZMfUc7fDmw/cDpS4AVU1zzJOAkgCRXDfTN1TziZzd/+dnNb35+85ef3fw1KWAOparG0Q912qzVstsrAJJsChxcVfcmWQ7sPencS9s1l0wqf1x+G7exdTlIskmSzSbWgZcA1wNLgYkn4w4DPtPWlwKHttEO9gLuneiaIEmSpDV2JbBLkp2TbAi8li5/PSLJ4iQT+fBY4JS2fgHwkiRbtIfBXgJc0LLaj5Ls1UY3OJRHs91aM84W2qcC57WRG9YHPllV/57kSuCcJEcA3wVe1Y4/H9gfWAbcD7xujHWTJElaUKrqoSRvoAuni4BTquqGJMcDV1XVUrpW2HcnKbouB69v596V5G/pQjHA8VV1V1v/M+BU4AnA59qyVqXm8bRqSY4c6BeiecTPbv7ys5vf/PzmLz+7+cvPbvzmdaCVJEmS5mIcWkmSJGnWzNtAO93UbZpbSbZPckmSm5LckOSvWrlTH88DSRYl+WqSz7btnZNc0T63s9vDBCTZqG0va/t3mst6C5JsnuTcJF9v37/f8Hs3PyT56/bfy+uTnJlkY797/ZTklCR3JLl+oGzk71lWMZWsRjcvA22Gm7pNc+sh4M1V9UxgL+D17TNy6uP54a+Amwa23wuc2D63u4EjWvkRwN1V9QvAie04za33A/9eVc8Afo3uc/R713NJtgP+Etijqp5F98DOa/G711en8vjpXUf6nmX1U8lqRPMy0DLc1G2aQ1V1e1V9pa3/iO7/VLfDqY97L8kS4HeAj7btAPsA57ZDJn9uE5/nucCL2vGaA0meBDwfOBmgqn5WVffg926+WB94QpL1gScCt+N3r5eq6jLgrknFo37PppxKdvy1XzfN10DrNLnzSPtT2O7AFTj18XzwD8DfABNzKG4F3FNVD7Xtwc/mkc+t7b+3Ha+58XTgh8DHWpeRj6YbB9zvXc9V1feA99ENZ3k73XfpavzuzSejfs/8/s2i+RpoezHNmqaXbpaRTwFvrKr7VnfoFGV+pmtZkpcBd1TV1YPFUxxaQ+zT2rc+8Gzgw1W1O/ATHv2z51T8/Hqi/an5QGBnYFtgE7o/VU/md2/+6fWUseuK+Rpoh5omV3MryQZ0YfaMqvp0K/7BxJ80M4OpjzV2zwMOSPJtuq48+9C12G7e/gwKj/1sHvnc2v4n8/g/w2ntWQ4sr6or2va5dAHX713/vRi4tap+WFU/Bz4N/CZ+9+aTUb9nfv9m0XwNtNNO3aa51fpynQzcVFV/P7DLqY97rKqOraolVbUT3ffq4qr6feAS4JXtsMmf28Tn+cp2vC0Mc6Sqvg/cluSXWtGLgBvxezcffBfYK8kT238/Jz47v3vzx6jfsymnkl3blV5XzNuJFZLsT9dyNDF12wlzXCUNSPJbwP8DruPRvphvpetHew6wA23q4zadXoAP0nWIvx94XVVdtdYrrkck2Rs4qqpeluTpdC22WwJfBf6gqh5MsjHwcbo+0ncBr62qW+aqzoIku9E90LchcAvdNOLr4feu95K8E3gN3SgxXwX+mK5Ppd+9nklyJt0UsYuBH9CNVvD/MeL3LMkf0f1/I8AJVfWxtfk+1iXzNtBKkiRJMH+7HEiSJEmAgVaSJEnznIFWkiRJ85qBVpIkSfOagVaSJEnzmoFWkiRJ85qBVtKClmSnJD9Ncs2I570mybIknx1X3SRJwzHQShJ8q6p2G+WEqjqbbuB7SdIcM9BKWmcleU6Sa5NsnGSTJDckedY05+yU5OtJPprk+iRnJHlxkv9McnOSPddW/SVJw1l/risgSeNSVVcmWQq8C3gC8Imqun6IU38BeBVwJHAl8HvAbwEH0E1TedB4aixJmgkDraR13fF0ofQB4C+HPOfWqroOIMkNwEVVVUmuA3YaSy0lSTNmlwNJ67otgU2BzYCNhzznwYH1lQPbK7EhQJJ6x0AraV13EvA/gTOA985xXSRJY2BLg6R1VpJDgYeq6pNJFgFfTLJPVV0813WTJM2eVNVc10GS5kySnYDPVtVqRz9Yxbl7A0dV1ctmuVqSpBHY5UDSQvcw8OSZTKwA/BNw91hqJUkami20kiRJmtdsoZUkSdK8ZqCVJEnSvGaglSRJ0rxmoJUkSdK8ZqCVJEnSvPb/A3pOGATfWr7qAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Define P-wave velocity for homogeneous medium\n",
    "vp0 = 1.0 # Vp [m/s]\n",
    "vp = vp0*np.ones((NZ+1,NX+1))\n",
    "\n",
    "# Plot Vp-model\n",
    "cax = plt.pcolormesh(x, z, vp)\n",
    "plt.title(\"Homogeneous P-wave velocity model\" )\n",
    "plt.xlabel(\"x [m]\")\n",
    "plt.ylabel(\"z [m]\")\n",
    "cbar = plt.colorbar(cax, orientation='vertical', pad=0.02)\n",
    "cbar.set_label('Vp [m/s]', labelpad=10)\n",
    "plt.gca().invert_yaxis()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hmm, looks a little bit boring, so let's build a more interesting model - e.g. a layer over a half-space:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqgAAAFNCAYAAADM2Z3wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de5RsZX3u++/DRVFREVGDCwgk4iXxBDWIeEyUiFuBqJijJuKNGLIxR7OjObqj8WRrkq17xFzEOIziiiCoiHoAFRlGQ0AgbgUFIVyCkSVeWIIgchUEw1q/88d8G4umV3d1r67Vb63+fhw1uuac75z1VhWFD+9tpqqQJEmSerHNSldAkiRJGmVAlSRJUlcMqJIkSeqKAVWSJEldMaBKkiSpKwZUSZIkdcWAKm2lknwnybNWuh6r1XJ8/kn+Kcnhy1WnSUnyu0m+NGbZ45K8fdJ1kjTdDKjSFEhy3yTHJPlukluTXJjk4JWu1zhaIPlpkh8nuSHJ6Ukeu9L1mgZVdXBVHQ+LC4GSNO0MqNJ02A64CngG8GDgfwCfTLLnCtbpXpJst4lDf11VOwK7AdcBx22xSkmSpo4BVZoCVXVbVf15VX2nqjZW1WnAt4FfHef8JPsl+UqSm5Jck+S9Se7Tjv1Dkr+bVf6zSV7fnj8yyclJfpjk20n+aKTcnyc5KclHk9wC/O4C7+N24GPA4zdRz+OTvKE9X5OkkrymbT+qtcAmyUOSnNbqdGN7vlsr95Ik58+67h8nObU9v2+Sv03yvSTXJjk6yf3mqMt92+f1+JF9D0vykyQPb9vPTXJRK/flJL+yifd13yTvTnJ1e7w7yX1Hjh/arnNLkm8lOajtPyvJ7yd5HHA08NTWEn1Tkie3+m83cp0XJrloE3U4Lsn72rCBHyf530l+rtXlxiTfSPLEkfKPa69/U5LLkjx/5NhDk5za6vtV4BdnvdZjW0v5DUn+I8lvz1UnSdoUA6o0hZI8Ang0cNmYp2wA/hjYBXgqcCDwmnbseOCwJNu0a+/Sjp/Y9n0W+DdgTdv/+iTPGbn2ocBJwE7ACQvUe0fgZcCFmyhyNnBAe/4M4Mr2F+DpwL/WcH/mbYAPAT8P7AH8BHhvK3cq8Jgke49c96UMwRjgnQyf3ROAR7X39dbZFamqO4FTgMNGdv82cHZVXZfkScCxwKuBhwIfAE4dDZ4j/l9g//aa+wD7AX/WPpP9gA8D/53hM3w68J1Zdbkc+APgK1W1Y1XtVFVfA34E/JeRoi8HPjLH64/W/88Y/jm4E/gK8PW2fRLwrlan7Rm+938GHg78N+CEJI9p1/kH4A5gV+D32oN27gOA0xk+74e3z+99SX55nnpJ0j0YUKUp08LDCcDxVfWNcc6pqguq6tyququqvsMQpp7Rjn0VuJkhfAK8BDirqq4Fngw8rKr+sqp+WlVXAv/Yysz4SlV9urXs/mQTVXhjkpuAdcCObLql9Wzg11swfjrw18DT2rFntONU1Y+q6uSqur2qbgXeMfJ+bgc+QwuWLag+liE8BvivwB9X1Q3t3P816/2M+hj3DKijQfe/Ah+oqvOqakMbK3onQxCd7WXAX1bVdVX1Q+AvgFe0Y0cAx1bV6e0z/P643yvDf1y8vL3PnYHnjNRvLp9q/yzcAXwKuKOqPlxVG4BPADMtqPszfE9/1b73M4HTGP5DZlvghcBbW8v+pa0eM54LfKeqPtT+efs6cDLwojHfkyQZUKVp0oLbR4CfAn84sn+m2/bHSV42x3mPbt3gP2hd8f+LodVsxt1Bh3u2wv088MjWzXtTC5lvAR4xcu5VY1T9b1ur389V1fOr6lutXj8eeezR9v+YoaXx1xlC0dWt5e7ugJrk/kk+kGHS2C3AOcBOLTzBPYPlS4FPt+D6MOD+wAUj7+fzbf9czgTul+QpSX6+1etTI5/NG2Z9NrsDj5zjOo8Evjuy/d2RcrsD31r4I5zTR4HntZbp32ZoYb5mnvLXjjz/yRzbO47U96qq2jirzmsYPquZMdGjx2b8PPCUWZ/Ly4CfG/9tSVrtNjWhQVJnWuvfMQzh8JCq+s+ZY1W10Iz+9zN0qx9WVbdmGF862qL1UeDSJPsAjwM+3fZfBXy7qvZm02px72TkxGHi1Gxnt7rdp6q+n+Rs4JXAQ4CZ8ZVvAB4DPKWqfpDkCQzvL+34PwO7tP2HMQxvALieIYj9clV9f4z6bUzyyXaNa4HTWqsrDJ/NO6rqHWO81asZgtvMkIw92r6Z6/ziXCfNrs4c9ft+kq8Av8XQIvv+Ma4zjquB3ZNsMxJS9wC+CfwQuIshWH9j5NiMqxiGQYwOPZCkRbEFVZoe72cIj8+bpyt9Ux4I3AL8OMMST//36MGqWg98jaHl9OSR638VuCXJm5LcL8m2SR6f5Mmb9U7mdzZD6/A5bfsshjGQX2pd0TPv5yfATa1r+22z3s9dDGMq/wbYmWFMJC1s/SNw1MhEpzWzxtTO9jHgdxhaAUe7z/8R+IPWupokD0jym0keOMc1TgT+LMMkq10Yxrx+tB07BnhVkgOTbNPqM9cyXNcCu6VNbhvxYeBPgP+Dn7Xubq7zgNuAP0myfZIDgOcBH2/fwSnAn7eW7F8CRtdqPQ14dJJXtHO3bxO6HrdMdZO0ChhQpSnQupdfzdDF/IP5uvM34Y0MXd23MgSrT8xR5niGkHP3JJsWRp7XXvfbDC2QH2RY6mpSzmYIoDMB9UsM3fLnjJR5N3C/Vp9zGbrpZ/sY8Czg/2uBdcabGMbCntuGB/wLQ2vsnKpqJqw9Evinkf3nM4xDfS9wY7vm727iMm8HzgcuBi5hmJj09nadrwKvAo5iGAt8NkNr62xnMrTA/iDJ9SP7P9XKf6qqbtvU+1iMqvop8HzgYIbP+H3AK0fGxv4hw3CAHzAsGfahkXNvBZ7NMK736lbmncBck8ckaU4ZJsRKWu2SPJ2hVW/PWWMP1bkk3wJeXVX/stJ1kaTlYAuqpJmVAV4HfNBwOl2SvJBhfOqZK10XSVouXQXUJAe1RZ3XJXnzStdHWg3a2MCbGNa0fPcKV0eLkOQshrHJr/U/LCSNI8nuSb6Y5PIMN+F43RxlkuQ9LY9d3NZ93rL17KWLvy0P802GRadnJmwcVlX/vqIVkyRJ2kok2RXYtaq+3iZ1XgC8YDRvJTmEYXLqIcBTgL+vqqdsyXr21IK6H7Cuqq5sA/Q/znCHGkmSJC2Dqrqm3UBjZlLj5QxrHI86FPhwDc5lWGd61y1Zz54C6hruufDzeu79gUmSJGkZJNmT4Q5y5806tOKZrKeF+jPHvnuNP0hyJHAkwLZs+6v350GTrpckSdKS3cqN11fVpu5YN6fn/MYD6kc3bFi44CwXXHznZcAdI7vWVtXa2eXaHehOBl5fVbfMPjzHpbfomNCeAup6hjuTzNiNn91p5W7tQ14L8KDsXE/JgbOLSJIkdeNf6qTvLlzqnn50wwa++oU9Fi44y7a7XnFHVe07X5m2csvJwAlVdcocRcbKZJPUUxf/14C9k+zV7pTyEuDUFa6TJEnSFlfAxiX8byEjt82+vKretYlipwKvbLP59wdurqprlu3NjaGbFtSquivJHwJfALYFjq2qyxY4TZIkaStUbJjM6nFPA14BXJLkorbvLcAeAFV1NPA5hhn864DbGe52t0V1E1ABqupzDB+KJEnSqjW0oC7/sM+q+hJzjzEdLVPAa5f9xRehq4AqSZKkwThd9lsrA6okSVJnimJDJzdTWgkGVEmSpA5Noot/WhhQJUmSOlPABgOqJEmSemILqiRJkrpR4BhUSZIk9WX1zuE3oEqSJHWnKMegSpIkqSMFG1ZvPjWgSpIk9Wa4k9TqZUCVJEnqTtgw/x1Jt2oGVEmSpM4UsNEufkmSJPVkNbegbrPSFZAkSZJG2YIqSZLUmeFWp6u3BdWAKkmS1KGNZUCVJElSJ2xBlSRJUleKsGEVTxUyoEqSJHXILn5JkiR1wy5+SZIkdSZsKLv4JUmS1IkCNjoGVZIkST2xi1+SJEndqLKLX5IkSZ3ZaAuqJEmSejHM4rcFVZIkSd2wi1+SJEkdcRa/JEmSurPBO0lJkiSpF0VW9RjU1fvOJUmS1CVbUCVJkjq00UlSkiRJ6oXLTEmSJKkrRZwkJUmSpL64zJQkSZK6UYUL9UuSJKknYSN28UuSJKkThS2okiRJ6oyz+CVJktSNImx0Fr8kSZJ6YguqJEmSulF4JylJkiR1JWxwFr8kSZJ6YQuqJEmSumMLqiRJkrpRlVXdgjqxd57k2CTXJbl0ZN/OSU5PckX7+5C2P0nek2RdkouTPGlS9ZIkSZoGG2qbRT8WMlc+m3X8wUk+m+TfklyW5FXL/sbGMMlofhxw0Kx9bwbOqKq9gTPaNsDBwN7tcSTw/gnWS5IkabU6jnvns1GvBf69qvYBDgD+Lsl9tkC97mFiAbWqzgFumLX7UOD49vx44AUj+z9cg3OBnZLsOqm6SZIk9ayAjWTRjwWvO3c+m/3SD0wSYMdW9q7leE+LsaXHoD6iqq4BqKprkjy87V8DXDVSbn3bd83sCyQ5kqGVlR24/2RrK0mStCIyVpf9HHZJcv7I9tqqWruI898LnApcDTwQ+J2q2riUimyOXiZJzRX5a66C7UNeC/Cg7DxnGUmSpGk2LDO1pFn811fVvpvx0s8BLgKeCfwicHqSf62qWzbjmou2paeHXTvTdd/+Xtf2rwd2Hym3G0NylyRJWpU2sM2iH8vgVcApbdjlOuDbwGOX48KLsaUD6qnA4e354cBnRva/ss3m3x+4eWYogCRJ0mpThI21+Mcy+B5wIECSRwCPAa5cjgsvxsS6+JOcyDD7a5ck64G3AX8FfDLJEQwfwItb8c8BhwDrgNsZ0rskSdKqtXEC7YibyGfbA1TV0cD/BI5LcgnDEMw3VdX1y16RBUwsoFbVYZs4dOAcZYthWQNJkqRVrwo2LE+L6KzrbjKfzRy/Gnj2sr/wIvUySUqSJEkjlqnLfioZUCVJkjozjEFdvbc6NaBKkiR1aMMYC+9vrQyokiRJndmMdVC3CgZUSZKk7tjFL0mSpM5stItfkiRJvZjUMlPTwoAqSZLUIbv4JUmS1I2ZW52uVqs3mkuSJKlLtqBKkiR1yElSkiRJ6obroEqSJKk7TpKSJElSP2p1T5IyoEqSJHWmcAyqJEmSOmMLqiRJkrrhJClJkiR1x4AqSZKkbqz2O0kZUCVJkjrkJClJkiT1o+zilyRJUkecJCVJkqTuGFAlSZLUDSdJSZIkqTtlQJUkSVJPnMUvSZKkbtQqn8W/zUpXQJIkSRplC6okSVKHHIMqSZKkjjiLX5IkSZ2xBVWSJEnd8E5SkiRJ6ksNM/lXKwOqJElSh1wHVZIkSd0oHIMqSZKkrjiLX5IkSZ1xDKokSZK6Yhe/JEmSulFlQJUkSVJnHIMqSZKkrjgGVZIkSV2xi1+SJEndKGJAlSRJUl9WcQ8/26x0BSRJkqRREwuoSXZP8sUklye5LMnr2v6dk5ye5Ir29yFtf5K8J8m6JBcnedKk6iZJktS1tszUYh8LSXJskuuSXDpPmQOSXNTy29nL+r7GNMkW1LuAN1TV44D9gdcm+SXgzcAZVbU3cEbbBjgY2Ls9jgTeP8G6SZIk9a2W8FjYccBBmzqYZCfgfcDzq+qXgRcvtfqbY2IBtaquqaqvt+e3ApcDa4BDgeNbseOBF7TnhwIfrsG5wE5Jdp1U/SRJkno2iRbUqjoHuGGeIi8FTqmq77Xy1y3Pu1mcLTIGNcmewBOB84BHVNU1MIRY4OGt2BrgqpHT1rd9kiRJq85wN6nFPZbBo4GHJDkryQVJXrksV12kic/iT7IjcDLw+qq6Jdlkup/rwL0+6iRHMgwBYAfuv1zVlCRJ6kax5HVQd0ly/sj22qpau4jztwN+FTgQuB/wlSTnVtU3l1KZpZpoQE2yPUM4PaGqTmm7r02ya1Vd07rwZ5qO1wO7j5y+G3D17Gu2D3ktwIOy82pegUGSJG2tClhaQL2+qvbdjFde365xG3BbknOAfYAtGlAnOYs/wDHA5VX1rpFDpwKHt+eHA58Z2f/KNpt/f+DmmaEAkiRJq80KdfF/Bvj1JNsluT/wFIZ5RFvUJFtQnwa8ArgkyUVt31uAvwI+meQI4Hv8bHbY54BDgHXA7cCrJlg3SZKkvk2gnzjJicABDEMB1gNvA7YHqKqjq+ryJJ8HLgY2Ah+sqk0uSTUpEwuoVfUl5h5XCsO4htnlC3jtpOojSZI0PSZzq9OqOmyMMn8D/M2yv/gieKtTSZKkHq3imTYGVEmSpN7UkmfxbxUMqJIkST2awhbUJO8Zo9gtVfVn8xUwoEqSJHVpKltQDwXeukCZNwMGVEmSpKkzhS2owFFVdfx8BZI8ZKGLGFAlSZJ6NIUBtarevRxlJrZQvyRJkpZo5k5Si310IslfJ3lQku2TnJHk+iQvH/d8A6okSZKW27Or6hbguQy3T3008N/HPdkufkmSpA4t061LV8r27e8hwIlVdUMyfguvAVWSJKlH0x1QP5vkG8BPgNckeRhwx7gn28UvSZLUoykcg5pkV4CqejPwVGDfqvpP4HaGJajGMm8LapJbFqoHcE1VPXrcF5QkSdLCMp0tqMe2ZaTOAj4PfAmgqm4Dbhv3Igt18X+rqp44X4EkF477YpIkSRpDMZVd/FV1cJIdgAOA3wL+Nsn3GMLq56vqe+NcZ6GA+sIxrjFOGUmSJI2tjy77paiqO2iBFCDJXsDBwHuT/FxV7bfQNeYNqFV15eh2kgeNnlNVN8wuI0mSpGUwhS2os7XseDPw8fb48TjnjTWLP8mrgb9kmIk183EV8AuLrqkkSZIWNsUBdVPZsarGyo7jLjP1RuCXq+r6xVdRkiRJizbFAZXNzI7jBtRvMSwPIEmSpEmbudXp9Nqs7DhuQP1T4MtJzgPunNlZVX+01BeWJEnSpk3pMlMzNis7jhtQPwCcCVwCbFxsDSVJkrRI0x1QNys7jhtQ76qq/2exF5ckSdKqtFnZcdyA+sUkRwKf5Z7NtDcs9YWXw6N/5Xa+8IV/W8kqSJIkzWvbXZd23pR38W9Wdhw3oL60/f3TkX0uMyVJkjQp0z1JarOy41gBtar2WmSlJEmStFRTeqvTGZubHbeZ72CSJy10gXHKSJIkaeu3XNlxoRbUDyU5AJivjfkY4IkLvZAkSZIWYTpbUJclOy4UUB8MXLDAi/xwgWtIkiRpkaZ0ktSyZMd5A2pV7bm4OkmSJGlZTGFAXa7sOO4sfkmSJG1JUxhQl4sBVZIkqTOpqe3iXxYGVEmSpB5N9zqom2XeZaZmJDkjySGz9q2dTJUkSZJ091qoi3l0JMn/leRdSf4uyW8t5tyxAiqwF/CmJG8b2bfvYl5IkiRJ45vp5l/MoxdJ3gf8AXAJcCnw6iT/MO7543bx3wQcCLwnyWeBly+2opIkSVqEjgLnEjwDeHxVFUCS4xnC6ljGbUFNVd1VVa8BTga+BDx8sTWVJEnSGJbQetpTCyrwH8AeI9u7AxePe/K4LahHzzypquOSXAK8dtwXkSRJ0iL1FTgX66HA5Um+2rafDHwlyakAVfX8+U4eK6BW1QdmbV8A/N7i6ypJkqSxTHdAfevmnOwyU5IkSR3qrMt+LEneC3ysqs7enOuMOwZVkiRJWsgVwN8l+U6SdyZ5wlIuYkCVJEnq0RSug1pVf19VT2WYxX8D8KEklyd5a5JHj3sdA6okSVJvpnwWf1V9t6reWVVPBF4K/BZw+bjnG1AlSZK0rJJsn+R5SU4A/gn4JvDCcc93kpQkSVKPOmoRHVeS/wIcBvwm8FXg48CRVXXbYq5jQJUkSerRFAZU4C3Ax4A3VtUNS72IAVWSJKkzoa8xpeOqqt9YjutMbAxqkh2SfDXJvyW5LMlftP17JTkvyRVJPpHkPm3/fdv2unZ8z0nVTZIkqXsTmMWf5Ngk1yW5dIFyT06yIcmLNuctLNUkJ0ndCTyzqvYBngAclGR/4J3AUVW1N3AjcEQrfwRwY1U9CjiqlZMkSVp9JjeL/zjgoPkKJNmWIYd9YbPfxxJNLKDW4Mdtc/v2KOCZwElt//HAC9rzQ9s27fiBSTKp+kmSJHVtAi2oVXUOw/qk8/lvwMnAdUur+Oab6DJTSbZNchHDGzwd+BZwU1Xd1YqsB9a052uAqwDa8ZuBh06yfpIkSd1agYX6k6xhWLP06M2/2tJNdJJUVW0AnpBkJ+BTwOPmKtb+ztVaeq+POsmRwJEAe6xxjpckSdo6LXGS1C5Jzh/ZXltVaxdx/ruBN1XVhpXsyN4iCa+qbkpyFrA/sFOS7Vor6W7A1a3YemB3YH2S7YAHM0cTdPuQ1wLsu88OUzi/TZIkaQxLSznXV9W+m/Gq+wIfb+F0F+CQJHdV1ac345qLNslZ/A9rLackuR/wLIZbXH0RmJkRdjjwmfb81LZNO35mVRlAJUnS6rOU7v1lSE1VtVdV7VlVezLMCXrNlg6nMNkW1F2B49tMsG2AT1bVaUn+nSGZvx24EDimlT8G+EiSdQwtpy+ZYN0kSZK6Nol1UJOcCBzAMBRgPfA2honsVNWKjjsdNbGAWlUXA0+cY/+VwH5z7L8DePGk6iNJkjRVJhBQq+qwRZT93eWvwXicZSRJktShabyT1HIxoEqSJPXIgCpJkqRuLNOkp2llQJUkSepMmHuB+NXCgCpJktSjVdyCOtFbnUqSJEmLZQuqJElSh5zFL0mSpL4YUCVJktQVA6okSZK6UXbxS5IkqTcGVEmSJPXEFlRJkiT1xYAqSZKkntiCKkmSpH4UtqBKkiSpMwZUSZIk9SLYxS9JkqTeGFAlSZLUk9TqTagGVEmSpN44SUqSJEm9cQyqJEmS+rKKA+o2K10BSZIkaZQtqJIkSR2yi1+SJEl9MaBKkiSpG2ULqiRJknpjQJUkSVIvvNWpJEmS+uOdpCRJktQTW1AlSZLUD291KkmSpN5k40rXYOUYUCVJknpkC6okSZJ64hhUSZIk9aNwFr8kSZL6YguqJEmS+mJAlSRJUi+8k5QkSZL6UrWqx6Bus9IVkCRJkkbZgipJktQhu/glSZLUFwOqJEmSemILqiRJkvpRwMbVm1ANqJIkST1avfnUgCpJktSj1dzFP/FlppJsm+TCJKe17b2SnJfkiiSfSHKftv++bXtdO77npOsmSZLUrZm1UBfzWECSY5Ncl+TSTRx/WZKL2+PLSfZZ9vc1hi2xDurrgMtHtt8JHFVVewM3Ake0/UcAN1bVo4CjWjlJkqRVKbX4xxiOAw6a5/i3gWdU1a8A/xNYu9lvZAkmGlCT7Ab8JvDBth3gmcBJrcjxwAva80PbNu34ga28JEnS6lJLfCx02apzgBvmOf7lqrqxbZ4L7LbUt7A5Jt2C+m7gT4CNbfuhwE1VdVfbXg+sac/XAFcBtOM3t/L3kOTIJOcnOf+HP9owybpLkiStiACpWvQD2GUmJ7XHkZtRjSOAf1qWN7RIE5skleS5wHVVdUGSA2Z2z1G0xjj2sx1Va2nNzfvus8MqHj4sSZK2ahsXLjKH66tq38196SS/wRBQf21zr7UUk5zF/zTg+UkOAXYAHsTQorpTku1aK+luwNWt/Hpgd2B9ku2ABzNPE7QkSdLWLGNMeprI6ya/wjA88+Cq+tFK1GFiXfxV9adVtVtV7Qm8BDizql4GfBF4USt2OPCZ9vzUtk07fmbVCn0zkiRJK2lCY1AXkmQP4BTgFVX1zc2/4tKsxDqobwI+nuTtwIXAMW3/McBHkqxjaDl9yQrUTZIkqQPjLRu1WElOBA5gGKu6HngbsD1AVR0NvJVhDtD72lz1u5ZjyMBibZGAWlVnAWe151cC+81R5g7gxVuiPpIkSb2bxEL9VXXYAsd/H/j95X/lxfFOUpIkST1axSMdt8RC/ZIkSdLYbEGVJEnqTUGWtszUVsGAKkmS1KNV3MVvQJUkSerR6s2nBlRJkqQerdRC/T0woEqSJPXIgCpJkqRuFOAkKUmSJPUilF38kiRJ6owBVZIkSV0xoEqSJKkbjkGVJElSbxyDKkmSpL4YUCVJktSPMqBKkiSpI4UBVZIkSZ1xkpQkSZJ6sponSW2z0hWQJEmSRtmCKkmS1KNV3IJqQJUkSepNARsNqJIkSeqGy0xJkiSpNwZUSZIkdcWAKkmSpG44BlWSJEl9KajVu1K/AVWSJKlHdvFLkiSpG3bxS5IkqTu2oEqSJKkrBlRJkiT1w4X6JUmS1JMCNjqLX5IkST2xBVWSJEldMaBKkiSpH+UyU5IkSepIQa3iO0lts9IVkCRJkkbZgipJktQju/glSZLUFSdJSZIkqRtVroMqSZKkztiCKkmSpJ6ULaiSJEnqR9mCKkmSpI4UzuKXJElSZ1yofzKSfCfJJUkuSnJ+27dzktOTXNH+PqTtT5L3JFmX5OIkT5pk3SRJknpVQG2sRT8WkuTYJNcluXQTx7vIY1viTlK/UVVPqKp92/abgTOqam/gjLYNcDCwd3scCbx/C9RNkiSpP1VDC+piHws7DjhonuNd5LGVuNXpocDx7fnxwAtG9n+4BucCOyXZdQXqJ0mStOIm0YJaVecAN8xTpIs8NumAWsA/J7kgyZFt3yOq6hqA9vfhbf8a4KqRc9e3fZIkSavPZFpQF9JFHpv0JKmnVdXVSR4OnJ7kG/OUzRz77vWfAi3ozoTdO7fd9Yo5x1Coe7sA1690JbQkfnfTze9vevndTa/HLPaEW7nxC/9SJ+2yhNfaYWbeT7O2qtYu4vyx8tikTTSgVtXV7e91ST4F7Adcm2TXqrqmNRlf14qvB3YfOX034Oo5rrkWWAuQ5PyRsa2aIn5308vvbrr5/U0vv7vpNSswjqWq5hsnOklj5bFJm1gXf5IHJHngzHPg2cClwKnA4a3Y4cBn2vNTgVe22WP7AzfPDAWQJEnSFtFFHptkC+ojgE8lmXmdj1XV55N8DfhkkiOA7wEvbuU/BxwCrANuB141wbpJkiStOklOBA4AdkmyHngbsD1AVR1NJ3lsYgG1qq4E9plj/4+AA+fYX8BrF/kyixlTob743U0vv7vp5vc3vfzuplc3311VHXGnRxIAAAVDSURBVLbA8aXksWWXWsX3eZUkSVJ/VmIdVEmSJGmTpjagJjkoyX+0W3G9eeEztCUl2T3JF5NcnuSyJK9r+73V7RRIsm2SC5Oc1rb3SnJe+94+keQ+bf992/a6dnzPlay3IMlOSU5K8o32+3uqv7vpkOSP278vL01yYpId/O31aa7bhS7ld5bk8Fb+iiSHz/Vaq9VUBtQk2wL/wHA7rl8CDkvySytbK81yF/CGqnocsD/w2vYdeavb6fA64PKR7XcCR7Xv7UbgiLb/CODGqnoUcFQrp5X198Dnq+qxDPMALsffXfeSrAH+CNi3qh4PbAu8BH97vTqOe98udFG/syQ7M0xQegrDMpxvmwm1mtKAyvBFrquqK6vqp8DHGW7NpU5U1TVV9fX2/FaG/5Ncg7e67V6S3YDfBD7YtgM8EzipFZn9vc18nycBB7byWgFJHgQ8HTgGoKp+WlU34e9uWmwH3C/JdsD9gWvwt9elTdwudLG/s+cAp1fVDVV1I3A69w69q9a0BtQubsOl8bSupycC5+GtbqfBu4E/AWbumfdQ4Kaquqttj343d39v7fjNrbxWxi8APwQ+1IZofLCtQ+3vrnNV9X3gbxmWX7yG4bd0Af72pslif2f+/uYxrQG1i9twaWFJdgROBl5fVbfMV3SOfX6nW1iS5wLXVdUFo7vnKFpjHNOWtx3wJOD9VfVE4DZ+1s04F7+/TrSu3UOBvYBHAg9g6Bqezd/e9NnUd+V3OI9pDahd3IZL80uyPUM4PaGqTmm7r53pQswSbnWriXsa8Pwk32EYOvNMhhbVnVq3I9zzu7n7e2vHH8y9u7205awH1lfVeW37JIbA6u+uf88Cvl1VP6yq/wROAf5P/O1Nk8X+zvz9zWNaA+rXgL3b7Mb7MAwkP3WF66QRbSzUMcDlVfWukUPe6rZjVfWnVbVbVe3J8Ls6s6peBnwReFErNvt7m/k+X9TK2wKwQqrqB8BVSR7Tdh0I/Dv+7qbB94D9k9y//ftz5rvztzc9Fvs7+wLw7CQPaS3oz277xBQv1J/kEIaWnW2BY6vqHStcJY1I8mvAvwKX8LOxjG9hGIf6SWAP2q1uq+qG9i/k9zIMEL8deFVVnb/FK667JTkAeGNVPTfJLzC0qO4MXAi8vKruTLID8BGGMcY3AC9pd5HTCknyBIYJbvcBrmS4TeE2+LvrXpK/AH6HYRWUC4HfZxiT6G+vMxm5XShwLcNs/E+zyN9Zkt9j+P9GgHdU1Ye25Pvo2dQGVEmSJG2dprWLX5IkSVspA6okSZK6YkCVJElSVwyokiRJ6ooBVZIkSV0xoEqSJKkrBlRJq1qSPZP8JMlFizzvd5KsS3LapOomSauVAVWS4FtV9YTFnFBVn2BYSF2StMwMqJK2WkmenOTiJDskeUCSy5I8foFz9kzyjSQfTHJpkhOSPCvJ/05yRZL9tlT9JWm12m6lKyBJk1JVX0tyKvB24H7AR6vq0jFOfRTwYuBI4GvAS4FfA57PcFvCF0ymxpIkMKBK2vr9JUPIvAP4ozHP+XZVXQKQ5DLgjKqqJJcAe06klpKku9nFL2lrtzOwI/BAYIcxz7lz5PnGke2N+B/2kjRxBlRJW7u1wP8ATgDeucJ1kSSNwZYASVutJK8E7qqqjyXZFvhykmdW1ZkrXTdJ0qalqla6DpK0YpLsCZxWVfPO7t/EuQcAb6yq5y5ztSRpVbOLX9JqtwF48FIW6gfeB9w4kVpJ0ipmC6okSZK6YguqJEmSumJAlSRJUlcMqJIkSeqKAVWSJEldMaBKkiSpK/8/buNjvZgGGg0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Define P-wave velocity for homogeneous medium\n",
    "vp1 = 2.0 # Vp in half-space [m/s]\n",
    "\n",
    "# Add half-space to model\n",
    "vp_2layer = vp\n",
    "\n",
    "vp_2layer[int((NZ+1)/2):NZ+1,:] = vp1\n",
    "\n",
    "# Plot Vp-model\n",
    "cax = plt.pcolormesh(x, z, vp_2layer)\n",
    "plt.title(\"2-layer P-wave velocity model\" )\n",
    "plt.xlabel(\"x [m]\")\n",
    "plt.ylabel(\"z [m]\")\n",
    "cbar = plt.colorbar(cax, orientation='vertical', pad=0.02)\n",
    "cbar.set_label('Vp [m/s]', labelpad=10)\n",
    "plt.gca().invert_yaxis()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let 's try an even more complex layer interface:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define 2-layer P-wave velocity model with sine-shape layer interface\n",
    "vp1 = 2.0                   # Vp in half-space [m/s]\n",
    "depth_layer = 250.0         # average depth of sine layer interface\n",
    "amp = 100.0                 # amplitude of sine layer undulation\n",
    "wave = 2.0 * np.pi / 100.0  # angular wavenumber of sine undulation\n",
    "\n",
    "# define function to create model with undulating sine-shaped interface\n",
    "def create_sin_model(NXtmp,NZtmp,DH):\n",
    "\n",
    "    # define homogeneous medium\n",
    "    model = vp0*np.ones((NZtmp+1,NXtmp+1))\n",
    "\n",
    "    # Add a layer interface undulating with a sine function to model\n",
    "    # loop over Cartesian grid\n",
    "    for j in range(NZtmp+1):\n",
    "        for i in range(NXtmp+1):\n",
    "        \n",
    "            depth = j*DH\n",
    "            x = i*DH\n",
    "        \n",
    "            # estimate depth of sine layer interface @ x\n",
    "            depth_sine = depth_layer + amp * np.sin(wave*x)\n",
    "        \n",
    "            if(depth>depth_sine):\n",
    "                model[j,i] = vp1 \n",
    "                \n",
    "    return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqgAAAFNCAYAAADM2Z3wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de5wlZX3n8c+XAURFQEQNDCAkYhLjChhUsmpCxCjgBWJkFTUSg8Gs7qobTbxsVjaJJnE38RYTyUQQNAi64AWNigQVNAgKglxEZcQLIxMQh5soGmZ++0dVy6Hpyzk9fbqfM/158zqvPlX1VNVTp/oMv/49l0pVIUmSJLViq+WugCRJkjTIAFWSJElNMUCVJElSUwxQJUmS1BQDVEmSJDXFAFWSJElNMUCVGpfkeUk+NYbj7pWkkmy92MfW/Bbr80/ywyQ/v1j1GpckJyV5w5Blv53kSeOuk6R2GaBKDUjy+CTnJ7klyYYk/5bk0QBVdUpVPXm56zhOfUDy4z7Yuj7Ju5Nsv9z1mgRVtX1VXQOjBYGS1DIDVGmZJdkB+Bjwd8DOwGrgz4CfLGe9xiGd2f7deXpVbQ88Cng08KdLVzNJUksMUKXl9zCAqjq1qjZW1Y+r6lNVdRlAkt9L8vmpwn2z8B8muTrJTUn+PkkGtv9+kqv6bWclecgwlUjywn6/25Jck+TFA9uuSPL0geVtktyYZL9++cA+A3xzkq8kOWig7GeTvDHJvwE/AuZsjq6q7wGfAB4xSz2/k+RX+/fP7z+Ph/fLL0ry4f79Y5J8oa/T+iTvSLJtv+34JH8z7bgfSfJH/fvdkpyR5PtJvpXkZbPU5cAk/55k1cC6304yde+2SvKaJN9M8oMkH0iy8yzH2i3JmX0GfW2SPxjYtirJ6/rj3Jbk4iR79NsqyUOTHAs8D/iTPhP90SR/nOSMaef5uyRvnaUO3+73uSzJ7UlOSPLgJJ/oz/uvSe4/UP4ZSa7sP+PPJvnlgW37J/lyv9/7ge2mnetpSS7t9z0/ySNnqpOklckAVVp+3wA2Jjk5yaGDAcAcnkaXZdwX+C/AUwCSHAG8Dngm8EDgc8CpQ9bjhv64OwAvBN6S5FH9tvcAzx8oexiwvqouTbIa+BfgDXQZ4FcBZyR54ED53wWOBe4HfGeuSvSB12HAJbMUORc4qH//68A1wG8MLJ/bv98I/A9gF+DXgIOBl/Tb3gc8eyqw7z/zJwOn9RnejwJfoctmHwy8IslTplekqi4AbgeeOLD6uf3xAV4GHNHXbzfgJuDvZ7muU4F1fblnAX+Z5OB+2x8BR9F9LjsAv08X7A/WZQ1wCvB/+mb/pwP/DBySZKf+OrcGng28d5Y6APwO8Ft0fzg9ne6PhdfRfY5b9ddEkof1dX4F3e/ax4GPJtm2/0Pgw/15dgb+X39c+n0fBZwIvBh4APCPwJlJ7jVHvSStIAao0jKrqluBxwMF/BPw/T6T9uA5dvvrqrq5qr4LfAbYr1//YuCvquqqqroT+Etgv2GyqFX1L1X1zeqcC3wKeEK/+Z+Bw9J1R4Au4JwKcp4PfLyqPl5Vm6rqbOAiumBqyklVdWVV3VlV/zFLFT6c5Gbg83RB5l/OUu5c7gpInwD81cDyb/TbqaqLq+qC/pzfpguCpsp9ju7znrq+ZwFfqKrr6AL/B1bVn1fVT/v+nf8EPGeW+pxKFzyS5H79dU/9UfBi4H9W1bqq+gnwv4FnZdrAqD4ofzzw6qq6o6ouBd5F9zkDvAj406r6en9/vlJVP5ilPj9TVeuB84Aj+1WHADdW1cVz7PZ3VXV9n8n+HHBhVV3S1/9DwP59uWcD/1JVZ/f39G+AewP/GTgQ2AZ4a1X9R1WdDnxp4Bx/APxjVV3YtxqcTNel5cD5rknSymCAKjWgDyh/r6p2p2va3g2YsRm29+8D738ETA0oegjwtr7Z9GZgAxC6TOCc+uztBX0T8810gdYuff2uA/4N+J0+G3coXbZu6pxHTp2z3/fxwK4Dh792vvMDR1TVTlX1kKp6SVX9OMkT+ubqHya5si93LvCEJD8HrALeDzwuyV7AjsCl/fU8LMnH+ib4W+kC3qnrKeA0+sCSLus5eD27Tbue1wGz/cHwPuCZffbvmcCXq+o7A8f60MBxrqLL7E4/1m7Ahqq6bWDdd7jrvu0BfHO+D3AWJ3NX9vv5zJ09Bbh+4P2PZ1ie+l3bjYFseFVtorvPq/tt3+s/5ymDmfOHAK+c9hnv0e8nSQaoUmuq6mvASczSB3Me1wIv7gO9qde9q+r8uXbqg6sz6LJgD66qneiabDNQbCrQOZIu2/i9gXO+d9o571tVfz14WQu4Fqrqc31z9fZV9Sv9urV0QfnLgPP6oO7f6boQfL4PlADeCXwN2KeqdqALMgev51S6bOZDgMf21z91Pd+adj33q6rBjPBgHb9KF3wdyt2b96eOdei0Y2038NlNuQ7Yuc/ATtkTGPyMf2HeD2zmz/nDwCOTPIKuC8cpM5RZiOvoAk2gGwBHF2R+D1gPrJ7qQtHbc+D9tcAbp30u96mqYbujSNrCGaBKyyzJLyV5ZZLd++U96DJ7FyzgcMcDr03yK/2xdkxy5Dz7AGwL3Av4PnBnkkPp+mQO+jDdCPuX0/VJnfLPwNOTPKUfzLNdkoOmrmdMzgX+G3f1N/3stGXo+rveCvwwyS8B/3XwAFV1Cd31vgs4q6pu7jd9Ebg1yauT3Lu/pkekn/ZrFu+jC5h/na6/5ZTjgTdOdbFI8sAkh0/fuaquBc4H/qr//B4JHMNdweS7gL9Isk86j0zygBnqcT3TBqFV1R3A6X0dv9h3C1kMHwCemuTgJNsAr6Rrpj8f+AJwJ/CyJFsneSbwmIF9/wn4wySP7a/nvkmeOi1Al7SCGaBKy+82ugzehUlupwtMr6D7H/5IqupDwJvoBvvc2h/n0CH2u40uwPoA3UCe5wJnTivzY7os497ABwfWXwscTpeh/D5dduyPGe+/L+fSBaDnzbIM3WCt59J9vv9E1xVgulOBJzGQ9ayqjXSDg/YDvgXcSBcg7jhHfU6lG7j16aq6cWD92+g+x08luY3u3j52lmMcBexFl5n8EHBc358X4M109+ZTdEH3CXT9Pac7AXh432z+4YH1JwP/ifmb94dWVV+ny6j/Hd1n9HS6qcJ+WlU/pevu8Ht0v0/P5u6/MxfR9UN9R799bV9WkgDI3bsISdLskrweeFhVPX/ewmpGkj3pujv8XD8oT5Ka5iMOJQ0l3fydx3DXyHJNgH7arD8CTjM4lTQpmmriT3JIkq+nm6T6NctdH0mddJPGXwt8oqrOm6+82pDkvnRdAn4LOG6ZqyOpAUn2SPKZdA9muTLJy2cokyRv7+Oxy3LXnNhLV89WmvjTPYnlG3T/kK6jmzPvqH6ErCRJkjZTkl2BXavqy/3AxIvppvn76kCZw4D/Tjfd4GOBt1XVbP3nx6KlDOpjgLVVdU3fwf40uoEXkiRJWgRVtb6qvty/v41ufubpc2UfDrynfzDIBcBOfWC7ZFoKUFdz98m81zHE5OKSJEkaXf+Ak/2BC6dtWvaYrKVBUplh3T36HyQ5lm5Cblax6lfvww732EmSJKkVt3HTjVX1wFH2ecpv3rd+sGHjyOe6+LKfXAncMbBqTVWtmV4uyfZ0Uwe+YoYBlEPFZOPUUoC6ju4pJFN2p5sP8G76D3kNwA7ZuR6bg5emdpIkSQvwr3X6d+YvdXc/2LCRL5615/wFp1m169V3VNUBc5XpH65xBnBKVX1whiJDxWTj1FIT/5eAfZLsnWRb4DlMmyhckiRpJShg0wL+m0//COITgKuq6s2zFDsTeEE/mv9A4JaqWr9oFzeEZjKoVXVnkv8GnAWsAk6sqiuXuVqSJEnLoNhY8wecC/A4uvmsL09yab/udcCeAFV1PPBxuhH8a4EfAS8cR0Xm0kyAClBVH6f7UCRJklasLoO6+N0+q+rzzNzHdLBMAS9d9JOPoKkAVZIkSZ1hmuy3VAaokiRJjSmKjY08TGk5GKBKkiQ1aBxN/JPCAFWSJKkxBWw0QJUkSVJLzKBKkiSpGQX2QZUkSVJbVu4YfgNUSZKk5hRlH1RJkiQ1pGDjyo1PDVAlSZJa0z1JauUyQJUkSWpO2Dj3E0m3aAaokiRJjSlgk038kiRJaslKzqButdwVkCRJkgaZQZUkSWpM96jTlZtBNUCVJElq0KYyQJUkSVIjzKBKkiSpKUXYuIKHChmgSpIkNcgmfkmSJDXDJn5JkiQ1Jmwsm/glSZLUiAI22QdVkiRJLbGJX5IkSc2osolfkiRJjdlkBlWSJEmt6Ebxm0GVJElSM2zilyRJUkMcxS9JkqTmbPRJUpIkSWpFkRXdB3XlXrkkSZKaZAZVkiSpQZscJCVJkqRWOM2UJEmSmlLEQVKSJElqi9NMSZIkqRlVOFG/JEmSWhI2YRO/JEmSGlGYQZUkSVJjHMUvSZKkZhRhk6P4JUmS1BIzqJIkSWpG4ZOkJEmS1JSw0VH8kiRJaoUZVEmSJDXHDKokSZKaUZUVnUEd25UnOTHJDUmuGFi3c5Kzk1zd/7x/vz5J3p5kbZLLkjxqXPWSJEmaBBtrq5Ff85kpPpu2fcckH03ylSRXJnnhol/YEMYZmp8EHDJt3WuAc6pqH+CcfhngUGCf/nUs8M4x1kuSJGmlOol7xmeDXgp8tar2BQ4C/jbJtktQr7sZW4BaVecBG6atPhw4uX9/MnDEwPr3VOcCYKcku46rbpIkSS0rYBMZ+TXvcWeOz6af+n5JAmzfl71zMa5pFEvdB/XBVbUeoKrWJ3lQv341cO1AuXX9uvXTD5DkWLosK9txn/HWVpIkaVlkqCb7GeyS5KKB5TVVtWaE/d8BnAlcB9wPeHZVbVpIRTZHK4OkZgr5a6aC/Ye8BmCH7DxjGUmSpEnWTTO1oFH8N1bVAZtx6qcAlwJPBH4BODvJ56rq1s045siWenjY9VNN9/3PG/r164A9BsrtThe5S5IkrUgb2Wrk1yJ4IfDBvtvlWuBbwC8txoFHsdQB6pnA0f37o4GPDKx/QT+a/0DglqmuAJIkSStNETbV6K9F8F3gYIAkDwZ+EbhmMQ48irE18Sc5lW701y5J1gHHAX8NfCDJMXQfwJF98Y8DhwFrgR/RRe+SJEkr1qYx5BFnic+2Aaiq44G/AE5KcjldF8xXV9WNi16ReYwtQK2qo2bZdPAMZYtuWgNJkqQVrwo2Lk5GdNpxZ43PprZfBzx50U88olYGSUmSJGnAIjXZTyQDVEmSpMZ0fVBX7qNODVAlSZIatHGIife3VAaokiRJjdmMeVC3CAaokiRJzbGJX5IkSY3ZZBO/JEmSWjGuaaYmhQGqJElSg2zilyRJUjOmHnW6Uq3c0FySJElNMoMqSZLUIAdJSZIkqRnOgypJkqTmOEhKkiRJ7aiVPUjKAFWSJKkxhX1QJUmS1BgzqJIkSWqGg6QkSZLUHANUSZIkNWOlP0nKAFWSJKlBDpKSJElSO8omfkmSJDXEQVKSJElqjgGqJEmSmuEgKUmSJDWnDFAlSZLUEkfxS5IkqRm1wkfxb7XcFZAkSZIGmUGVJElqkH1QJUmS1BBH8UuSJKkxZlAlSZLUDJ8kJUmSpLZUN5J/pTJAlSRJapDzoEqSJKkZhX1QJUmS1BRH8UuSJKkx9kGVJElSU2zilyRJUjOqDFAlSZLUGPugSpIkqSn2QZUkSVJTbOKXJElSM4oYoEqSJKktK7iFn62WuwKSJEnSoLEFqEn2SPKZJFcluTLJy/v1Oyc5O8nV/c/79+uT5O1J1ia5LMmjxlU3SZKkpvXTTI36mk+SE5PckOSKOcoclOTSPn47d1Gva0jjzKDeCbyyqn4ZOBB4aZKHA68BzqmqfYBz+mWAQ4F9+texwDvHWDdJkqS21QJe8zsJOGS2jUl2Av4BeEZV/Qpw5EKrvznGFqBW1fqq+nL//jbgKmA1cDhwcl/sZOCI/v3hwHuqcwGwU5Jdx1U/SZKklo0jg1pV5wEb5ijyXOCDVfXdvvwNi3M1o1mSPqhJ9gL2By4EHlxV66ELYoEH9cVWA9cO7LauXydJkrTidE+TGu21CB4G3D/JZ5NcnOQFi3LUEY19FH+S7YEzgFdU1a3JrNH9TBvu8VEnOZauCwDbcZ/FqqYkSVIzigXPg7pLkosGltdU1ZoR9t8a+FXgYODewBeSXFBV31hIZRZqrAFqkm3ogtNTquqD/errk+xaVev7Jvyp1PE6YI+B3XcHrpt+zP5DXgOwQ3ZeyTMwSJKkLVUBCwtQb6yqAzbjzOv6Y9wO3J7kPGBfYEkD1HGO4g9wAnBVVb15YNOZwNH9+6OBjwysf0E/mv9A4JaprgCSJEkrzTI18X8EeEKSrZPcB3gs3TiiJTXODOrjgN8FLk9yab/udcBfAx9IcgzwXe4aHfZx4DBgLfAj4IVjrJskSVLbxtBOnORU4CC6rgDrgOOAbQCq6viquirJJ4HLgE3Au6pq1impxmVsAWpVfZ6Z+5VC169hevkCXjqu+kiSJE2O8TzqtKqOGqLM/wX+76KffAQ+6lSSJKlFK3ikjQGqJElSa2rBo/i3CAaokiRJLZrADGqStw9R7Naq+tO5ChigSpIkNWkiM6iHA6+fp8xrAANUSZKkiTOBGVTgLVV18lwFktx/voMYoEqSJLVoAgPUqnrrYpQZ20T9kiRJWqCpJ0mN+mpEkv+TZIck2yQ5J8mNSZ4/7P4GqJIkSVpsT66qW4Gn0T0+9WHAHw+7s038kiRJDVqkR5cul236n4cBp1bVhmT4DK8BqiRJUosmO0D9aJKvAT8GXpLkgcAdw+5sE78kSVKLJrAPapJdAarqNcCvAQdU1X8AP6Kbgmooc2ZQk9w6Xz2A9VX1sGFPKEmSpPllMjOoJ/bTSH0W+CTweYCquh24fdiDzNfE/82q2n+uAkkuGfZkkiRJGkIxkU38VXVoku2Ag4DfBv4myXfpgtVPVtV3hznOfAHq7wxxjGHKSJIkaWhtNNkvRFXdQR+QAiTZGzgUeEeSn6uqx8x3jDkD1Kq6ZnA5yQ6D+1TVhullJEmStAgmMIM6XR873gKc1r9+OMx+Q43iT/Ji4M/pRmJNfVwF/PzINZUkSdL8JjhAnS12rKqhYsdhp5l6FfArVXXj6FWUJEnSyCY4QGUzY8dhA9Rv0k0PIEmSpHGbetTp5Nqs2HHYAPW1wPlJLgR+MrWyql620BNLkiRpdhM6zdSUzYodhw1Q/xH4NHA5sGnUGkqSJGlEkx2gblbsOGyAemdV/dGoB5ckSdKKtFmx47AB6meSHAt8lLunaTcs9MST5qzrvrKk53vKbvsu6fm2dEt9/7ZkS/27uZT3bku+tpVgKe/flv7/BK9vca3adWH7TXgT/2bFjsMGqM/tf752YJ3TTEmSJI3LZA+S2qzYcagAtar2HrFSkiRJWqgJfdTplM2NHbeaa2OSR813gGHKSJIkacu3WLHjfBnUdyc5CJgrx3wCsP98J5IkSdIIJjODuiix43wB6o7AxfOc5PvzHEOSJEkjmtBBUosSO84ZoFbVXqPVSZIkSYtiAgPUxYodhx3FL0mSpKU0gQHqYjFAlSRJakxqYpv4F4UBqiRJUosmex7UzTLnNFNTkpyT5LBp69aMp0qSJEn62Vyoo7wakuSZSd6c5G+T/PYo+w4VoAJ7A69OctzAugNGOZEkSZKGN9XMP8qrFUn+AfhD4HLgCuDFSf5+2P2HbeK/GTgYeHuSjwLPH7WikiRJGkFDAecC/AbwiKoqgCQn0wWrQxk2g5qqurOqXgKcAXweeNCoNZUkSdIQFpA9bSmDCnwd2HNgeQ/gsmF3HjaDevzUm6o6KcnlwEuHPYkkSZJG1FbAOaoHAFcl+WK//GjgC0nOBKiqZ8y181ABalX947Tli4HfH72ukiRJGspkB6iv35ydnWZKkiSpQY012Q8lyTuA91XVuZtznGH7oEqSJEnzuRr42yTfTvKmJPst5CAGqJIkSS2awHlQq+ptVfVrdKP4NwDvTnJVktcnediwxzFAlSRJas2Ej+Kvqu9U1Zuqan/gucBvA1cNu78BqiRJkhZVkm2SPD3JKcAngG8AvzPs/g6SkiRJalFDGdFhJfkt4CjgqcAXgdOAY6vq9lGOY4AqSZLUogkMUIHXAe8DXlVVGxZ6EANUSZKkxoS2+pQOq6p+czGOM7Y+qEm2S/LFJF9JcmWSP+vX753kwiRXJ3l/km379ffql9f22/caV90kSZKaN4ZR/ElOTHJDkivmKffoJBuTPGtzLmGhxjlI6ifAE6tqX2A/4JAkBwJvAt5SVfsANwHH9OWPAW6qqocCb+nLSZIkrTzjG8V/EnDIXAWSrKKLw87a7OtYoLEFqNX5Yb+4Tf8q4InA6f36k4Ej+veH98v02w9OknHVT5IkqWljyKBW1Xl085PO5b8DZwA3LKzim2+s00wlWZXkUroLPBv4JnBzVd3ZF1kHrO7frwauBei33wI8YJz1kyRJatYyTNSfZDXdnKXHb/7RFm6sg6SqaiOwX5KdgA8BvzxTsf7nTNnSe3zUSY4FjgXYc/XWnHXRVxaptm0567ot87o0+bbk380t+dpWgi35/m3J1wZb/vUt1AIHSe2S5KKB5TVVtWaE/d8KvLqqNi5nQ/aSjOKvqpuTfBY4ENgpydZ9lnR34Lq+2DpgD2Bdkq2BHZkhBd1/yGsADth3uwkc3yZJkjSEhUU5N1bVAZtx1gOA0/rgdBfgsCR3VtWHN+OYIxvnKP4H9plTktwbeBLdI64+A0yNCDsa+Ej//sx+mX77p6vKAFSSJK08C2neX4Soqar2rqq9qmovujFBL1nq4BTGm0HdFTi5Hwm2FfCBqvpYkq/SReZvAC4BTujLnwC8N8lauszpc8ZYN0mSpKaNYx7UJKcCB9F1BVgHHEc3kJ2qWtZ+p4PGFqBW1WXA/jOsvwZ4zAzr7wCOHFd9JEmSJsoYAtSqOmqEsr+3+DUYjk+SkiRJatAkPklqsRigSpIktcgAVZIkSc1YpEFPk8oAVZIkqTFh5gniVwoDVEmSpBat4AzqWB91KkmSJI3KDKokSVKDHMUvSZKkthigSpIkqSkGqJIkSWpG2cQvSZKk1higSpIkqSVmUCVJktQWA1RJkiS1xAyqJEmS2lGYQZUkSVJjDFAlSZLUimATvyRJklpjgCpJkqSWpFZuhGqAKkmS1BoHSUmSJKk19kGVJElSW1ZwgLrVcldAkiRJGmQGVZIkqUE28UuSJKktBqiSJElqRplBlSRJUmsMUCVJktQKH3UqSZKk9vgkKUmSJLXEDKokSZLa4aNOJUmS1JpsWu4aLB8DVEmSpBaZQZUkSVJL7IMqSZKkdhSO4pckSVJbzKBKkiSpLQaokiRJaoVPkpIkSVJbqlZ0H9StlrsCkiRJ0iAzqJIkSQ2yiV+SJEltMUCVJElSS8ygSpIkqR0FbFq5EaoBqiRJUotWbnxqgCpJktSildzEP/ZpppKsSnJJko/1y3snuTDJ1Unen2Tbfv29+uW1/fa9xl03SZKkZk3NhTrKax5JTkxyQ5IrZtn+vCSX9a/zk+y76Nc1hKWYB/XlwFUDy28C3lJV+wA3Acf0648BbqqqhwJv6ctJkiStSKnRX0M4CThkju3fAn6jqh4J/AWwZrMvZAHGGqAm2R14KvCufjnAE4HT+yInA0f07w/vl+m3H9yXlyRJWllqga/5Dlt1HrBhju3nV9VN/eIFwO4LvYTNMe4M6luBPwE29csPAG6uqjv75XXA6v79auBagH77LX35u0lybJKLklz0/R9sHGfdJUmSlkWAVI38AnaZipP617GbUY1jgE8sygWNaGyDpJI8Dbihqi5OctDU6hmK1hDb7lpRtYY+3XzAvtut4O7DkiRpi7Zp/iIzuLGqDtjcUyf5TboA9fGbe6yFGOco/scBz0hyGLAdsANdRnWnJFv3WdLdgev68uuAPYB1SbYGdmSOFLQkSdKWLEMMehrLeZNH0nXPPLSqfrAcdRhbE39Vvbaqdq+qvYDnAJ+uqucBnwGe1Rc7GvhI//7Mfpl++6erlunOSJIkLacx9UGdT5I9gQ8Cv1tV39j8Iy7McsyD+mrgtCRvAC4BTujXnwC8N8lauszpc5ahbpIkSQ0YbtqoUSU5FTiIrq/qOuA4YBuAqjoeeD3dGKB/6Meq37kYXQZGtSQBalV9Fvhs//4a4DEzlLkDOHIp6iNJktS6cUzUX1VHzbP9RcCLFv/Mo/FJUpIkSS1awT0dl2KifkmSJGloZlAlSZJaU5CFTTO1RTBAlSRJatEKbuI3QJUkSWrRyo1PDVAlSZJatFwT9bfAAFWSJKlFBqiSJElqRgEOkpIkSVIrQtnEL0mSpMYYoEqSJKkpBqiSJElqhn1QJUmS1Br7oEqSJKktBqiSJElqRxmgSpIkqSGFAaokSZIa4yApSZIktWQlD5LaarkrIEmSJA0ygypJktSiFZxBNUCVJElqTQGbDFAlSZLUDKeZkiRJUmsMUCVJktQUA1RJkiQ1wz6okiRJaktBrdyZ+g1QJUmSWmQTvyRJkpphE78kSZKaYwZVkiRJTTFAlSRJUjucqF+SJEktKWCTo/glSZLUEjOokiRJaooBqiRJktpRTjMlSZKkhhTUCn6S1FbLXQFJkiRpkBlUSZKkFtnEL0mSpKY4SEqSJEnNqHIeVEmSJDXGDKokSZJaUmZQJUmS1I4ygypJkqSGFI7ilyRJUmOcqH88knw7yeVJLk1yUb9u5yRnJ7m6/3n/fn2SvD3J2iSXJXnUOOsmSZLUqgJqU438mk+SE5PckOSKWbY3EY8txZOkfrOq9quqA/rl1wDnVNU+wDn9MsChwD7961jgnUtQN0mSpPZUdRnUUV/zOwk4ZI7tTcRjy/Go08OBk/v3JwNHDKx/T3UuAHZKsusy1E+SJGnZjSODWlXnARvmKNJEPDbuALWATyW5OMmx/boHV9V6gP7ng/r1q4FrB/Zd16+TJElaecaTQZ1PE/HYuAdJPa6qrkvyIODsJF+bo2xmWHePPwX6QHcq2P3Jql2vnrEPhZq3C3DjcldCC+K9m2zev8nlvQGMOQgAAAZOSURBVJtcvzjqDrdx01n/WqfvsoBzbTc17qe3pqrWjLD/UPHYuI01QK2q6/qfNyT5EPAY4Poku1bV+j5lfENffB2wx8DuuwPXzXDMNcAagCQXDfRt1QTx3k0u791k8/5NLu/d5JoWMA6lqubqJzpOQ8Vj4za2Jv4k901yv6n3wJOBK4AzgaP7YkcDH+nfnwm8oB89diBwy1RXAEmSJC2JJuKxcWZQHwx8KMnUed5XVZ9M8iXgA0mOAb4LHNmX/zhwGLAW+BHwwjHWTZIkacVJcipwELBLknXAccA2AFV1PI3EY2MLUKvqGmDfGdb/ADh4hvUFvHTE04zSp0Jt8d5NLu/dZPP+TS7v3eRq5t5V1VHzbF9IPLboUiv4Oa+SJElqz3LMgypJkiTNamID1CSHJPl6/yiu18y/h5ZSkj2SfCbJVUmuTPLyfr2Pup0ASVYluSTJx/rlvZNc2N+39yfZtl9/r355bb99r+WstyDJTklOT/K1/vv3a37vJkOS/9H/e3lFklOTbOd3r00zPS50Id+zJEf35a9OcvRM51qpJjJATbIK+Hu6x3E9HDgqycOXt1aa5k7glVX1y8CBwEv7e+SjbifDy4GrBpbfBLylv283Acf0648BbqqqhwJv6ctpeb0N+GRV/RLdOICr8HvXvCSrgZcBB1TVI4BVwHPwu9eqk7jn40JH+p4l2ZlugNJj6abhPG4qqNWEBqh0N3JtVV1TVT8FTqN7NJcaUVXrq+rL/fvb6P4nuRofddu8JLsDTwXe1S8HeCJwel9k+n2bup+nAwf35bUMkuwA/DpwAkBV/bSqbsbv3aTYGrh3kq2B+wDr8bvXpFkeFzrq9+wpwNlVtaGqbgLO5p5B74o1qQFqE4/h0nD6pqf9gQvxUbeT4K3AnwBTz8x7AHBzVd3ZLw/em5/dt377LX15LY+fB74PvLvvovGufh5qv3eNq6rvAX9DN/3ierrv0sX43Zsko37P/P7NYVID1CYew6X5JdkeOAN4RVXdOlfRGdZ5T5dYkqcBN1TVxYOrZyhaQ2zT0tsaeBTwzqraH7idu5oZZ+L9a0TftHs4sDewG3Bfuqbh6fzuTZ7Z7pX3cA6TGqA28RguzS3JNnTB6SlV9cF+9fVTTYhZwKNuNXaPA56R5Nt0XWeeSJdR3alvdoS735uf3bd++47cs9lLS2cdsK6qLuyXT6cLWP3ete9JwLeq6vtV9R/AB4H/jN+9STLq98zv3xwmNUD9ErBPP7pxW7qO5Gcuc500oO8LdQJwVVW9eWCTj7ptWFW9tqp2r6q96L5Xn66q5wGfAZ7VF5t+36bu57P68mYAlklV/TtwbZJf7FcdDHwVv3eT4LvAgUnu0//7OXXv/O5NjlG/Z2cBT05y/z6D/uR+nZjgifqTHEaX2VkFnFhVb1zmKmlAkscDnwMu566+jK+j64f6AWBP+kfdVtWG/h/kd9B1EP8R8MKqumjJK66fSXIQ8KqqelqSn6fLqO4MXAI8v6p+kmQ74L10fYw3AM/pnyKnZZJkP7oBbtsC19A9pnAr/N41L8mfAc+mmwXlEuBFdH0S/e41JgOPCwWupxuN/2FG/J4l+X26/zcCvLGq3r2U19GyiQ1QJUmStGWa1CZ+SZIkbaEMUCVJktQUA1RJkiQ1xQBVkiRJTTFAlSRJUlMMUCVJktQUA1RJK1qSvZL8OMmlI+737CRrk3xsXHWTpJXKAFWS4JtVtd8oO1TV++kmUpckLTIDVElbrCSPTnJZku2S3DfJlUkeMc8+eyX5WpJ3JbkiySlJnpTk35JcneQxS1V/SVqptl7uCkjSuFTVl5KcCbwBuDfwz1V1xRC7PhQ4EjgW+BLwXODxwDPoHkt4xHhqLEkCA1RJW74/pwsy7wBeNuQ+36qqywGSXAmcU1WV5HJgr7HUUpL0MzbxS9rS7QxsD9wP2G7IfX4y8H7TwPIm/MNeksbOAFXSlm4N8L+AU4A3LXNdJElDMBMgaYuV5AXAnVX1viSrgPOTPLGqPr3cdZMkzS5Vtdx1kKRlk2Qv4GNVNefo/ln2PQh4VVU9bZGrJUkrmk38kla6jcCOC5moH/gH4Kax1EqSVjAzqJIkSWqKGVRJkiQ1xQBVkiRJTTFAlSRJUlMMUCVJktQUA1RJkiQ15f8DfksD+tzU0P4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Create Vp-model\n",
    "vp_sine = create_sin_model(NX,NZ,dh)\n",
    "\n",
    "# Plot Vp-model\n",
    "cax = plt.pcolormesh(x, z, vp_sine)\n",
    "plt.title(\"Sine layer P-wave velocity model\" )\n",
    "plt.xlabel(\"x [m]\")\n",
    "plt.ylabel(\"z [m]\")\n",
    "cbar = plt.colorbar(cax, orientation='vertical', pad=0.02)\n",
    "cbar.set_label('Vp [m/s]', labelpad=10)\n",
    "plt.gca().invert_yaxis()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There seem to be something wrong with this model. While there is an undulating layer interface, it does not correctly follow the defined sine-shape. The reason is a spatial undersampling of the model. Let's try a finer model discretization:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NX =  200\n",
      "NZ =  100\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqgAAAFNCAYAAADM2Z3wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZxkdX3v/9ebARw3RETNsAmJGGP8CRgUctWEiJHFBaPyU5SIBC/mpzdq0MTl5koWTfTexC0mwYkgaBT1hxsaFXEDDQIysovKiAsjRERAcEHDzOf+cU5L0fRS1d3Vdarq9eRRj6lzzvec8/3W6Wo+/V1TVUiSJEldsdWoMyBJkiT1MkCVJElSpxigSpIkqVMMUCVJktQpBqiSJEnqFANUSZIkdYoBqtRxSZ6d5FNDuO7uSSrJ1it9bS1upT7/JD9O8qsrla9hSXJyktf0mfbbSR437DxJ6i4DVKkDkjw6yTlJfpTkhiT/keQRAFX17qp6/KjzOExtQPKzNtj6fpJ3JLnHqPM1DqrqHlV1FQwWBEpSlxmgSiOWZDvgY8A/AjsAOwN/Bfx8lPkahjTm+73zpKq6B/Bw4BHAX6xeziRJXWKAKo3egwCq6tSq2lxVP6uqT1XVJQBJnpvkizOJ22bhP05yZZIbk/xTkvQc/6MkV7THzkjygH4ykeTo9rxbklyV5Pk9xy5L8qSe7W2SXJ9k73Z7/7YG+KYkFyc5oCft55O8Nsl/AD8FFmyOrqrvAZ8AHjpPPr+T5Lfa90e2n8dD2u3nJflw+/6RSb7U5unaJG9Nsm177IQkfz/ruh9Jclz7fqckH0jygyTfSvKiefKyf5L/TLKmZ98fJJl5dlsleUWSbyb5YZL3J9lhnmvtlOT0tgZ9Y5L/3nNsTZJXtde5JcmGJLu2xyrJA5McCzwb+PO2JvqjSf4syQdm3ecfk7xpnjx8uz3nkiQ/SXJikvsn+UR7308nuXdP+icnubz9jD+f5Dd6ju2T5Cvtee8D1s661xOTXNSee06Sh82VJ0nTyQBVGr1vAJuTnJLkkN4AYAFPpKll3Av4f4GDAJI8BXgV8FTgvsAXgFP7zMd17XW3A44G3pjk4e2xdwJH9qQ9FLi2qi5KsjPw78BraGqAXwZ8IMl9e9L/IXAscE/gOwtlog28DgUunCfJWcAB7fvfAa4Cfrdn+6z2/WbgT4Edgd8GDgRe0B57D/CMmcC+/cwfD7y3reH9KHAxTW32gcBLkhw0OyNVdS7wE+CxPbuf1V4f4EXAU9r87QTcCPzTPOU6FdjUpns68LdJDmyPHQccQfO5bAf8EU2w35uX9cC7gf/dNvs/Cfg34OAk27fl3Bp4BvCuefIA8DTg92n+cHoSzR8Lr6L5HLdqy0SSB7V5fgnNz9rHgY8m2bb9Q+DD7X12AP7/9rq05z4cOAl4PnAf4G3A6UnuskC+JE0RA1RpxKrqZuDRQAH/CvygrUm7/wKnva6qbqqq7wKfA/Zu9z8f+LuquqKqbgP+Fti7n1rUqvr3qvpmNc4CPgU8pj38b8ChabojQBNwzgQ5RwIfr6qPV9WWqjoTuIAmmJpxclVdXlW3VdV/zZOFDye5CfgiTZD5t/OkO4vbA9LHAH/Xs/277XGqakNVndve89s0QdBMui/QfN4z5Xs68KWquoYm8L9vVf11Vf2i7d/5r8Az58nPqTTBI0nu2ZZ75o+C5wP/s6o2VdXPgb8Enp5ZA6PaoPzRwMur6taqugh4O83nDPA84C+q6uvt87m4qn44T35+qaquBc4GDm93HQxcX1UbFjjtH6vq+21N9heA86rqwjb/HwL2adM9A/j3qjqzfaZ/D9wV+G/A/sA2wJuq6r+q6jTgyz33+O/A26rqvLbV4BSaLi37L1YmSdPBAFXqgDagfG5V7ULTtL0TMGczbOs/e97/FJgZUPQA4M1ts+lNwA1AaGoCF9TW3p7bNjHfRBNo7djm7xrgP4CntbVxh9DU1s3c8/CZe7bnPhpY13P5qxe7P/CUqtq+qh5QVS+oqp8leUzbXP3jJJe36c4CHpPkV4A1wPuARyXZHbgXcFFbngcl+VjbBH8zTcA7U54C3ksbWNLUevaWZ6dZ5XkVMN8fDO8BntrW/j0V+EpVfafnWh/quc4VNDW7s6+1E3BDVd3Ss+873P7cdgW+udgHOI9TuL32+0gWrj0F+H7P+5/NsT3zs7YTPbXhVbWF5jnv3B77Xvs5z+itOX8A8NJZn/Gu7XmSZIAqdU1VfQ04mXn6YC7iauD5baA387prVZ2z0EltcPUBmlqw+1fV9jRNtulJNhPoHE5T2/i9nnu+a9Y9715Vr+st1hLKQlV9oW2uvkdV/Wa7byNNUP4i4Ow2qPtPmi4EX2wDJYB/Ab4G7FlV29EEmb3lOZWmNvMBwH5t+WfK861Z5blnVfXWCPfm8as0wdch3LF5f+Zah8y61tqez27GNcAObQ3sjN2A3s/41xb9wOb+nD8MPCzJQ2m6cLx7jjRLcQ1NoAk0A+BogszvAdcCO890oWjt1vP+auC1sz6Xu1VVv91RJE04A1RpxJI8OMlLk+zSbu9KU7N37hIudwLwyiS/2V7rXkkOX+QcgG2BuwA/AG5LcghNn8xeH6YZYf9imj6pM/4NeFKSg9rBPGuTHDBTniE5C/gf3N7f9POztqHp73oz8OMkDwb+v94LVNWFNOV9O3BGVd3UHjofuDnJy5PctS3TQ9NO+zWP99AEzL9D099yxgnAa2e6WCS5b5LDZp9cVVcD5wB/135+DwOO4fZg8u3A3yTZM42HJbnPHPn4PrMGoVXVrcBpbR7Pb7uFrIT3A09IcmCSbYCX0jTTnwN8CbgNeFGSrZM8FXhkz7n/Cvxxkv3a8tw9yRNmBeiSppgBqjR6t9DU4J2X5Cc0gellNP/DH0hVfQh4Pc1gn5vb6xzSx3m30ARY76cZyPMs4PRZaX5GU8u4B/DBnv1XA4fR1FD+gKZ27M8Y7u+Xs2gC0LPn2YZmsNazaD7ff6XpCjDbqcDj6Kn1rKrNNIOD9ga+BVxPEyDea4H8nEozcOuzVXV9z/4303yOn0pyC82z3W+eaxwB7E5TM/kh4Pi2Py/AG2iezadogu4Tafp7znYi8JC22fzDPftPAf4fFm/e71tVfZ2mRv0faT6jJ9FMFfaLqvoFTXeH59L8PD2DO/7MXEDTD/Wt7fGNbVpJAiB37CIkSfNL8mrgQVV15KKJ1RlJdqPp7vAr7aA8Seo0lziU1Jc083cew+0jyzUG2mmzjgPea3AqaVx0qok/ycFJvp5mkupXjDo/khppJo2/GvhEVZ29WHp1Q5K703QJ+H3g+BFnR1IHJNk1yefSLMxyeZIXz5EmSd7SxmOX5PY5sVcvn11p4k+zEss3aH6RbqKZM++IdoSsJEmSlinJOmBdVX2lHZi4gWaav6/2pDkU+BOa6Qb3A95cVfP1nx+KLtWgPhLYWFVXtR3s30sz8EKSJEkroKquraqvtO9voZmfefZc2YcB72wXBjkX2L4NbFdNlwLUnbnjZN6b6GNycUmSJA2uXeBkH+C8WYdGHpN1aZBU5th3p/4HSY6lmZCbNaz5rbux3Z1OkiRJ6opbuPH6qrrvIOcc9Ht3rx/esHnge2245OeXA7f27FpfVetnp0tyD5qpA18yxwDKvmKyYepSgLqJZhWSGbvQzAd4B+2HvB5gu+xQ++XA1cmdJEnSEny6TvvO4qnu6Ic3bOb8M3ZbPOEsa9ZdeWtV7btQmnZxjQ8A766qD86RpK+YbJi61MT/ZWDPJHsk2RZ4JrMmCpckSZoGBWxZwn+LaZcgPhG4oqreME+y04HntKP59wd+VFXXrljh+tCZGtSqui3J/wDOANYAJ1XV5SPOliRJ0ggUm2vxgHMJHkUzn/WlSS5q970K2A2gqk4APk4zgn8j8FPg6GFkZCGdCVABqurjNB+KJEnS1GpqUFe+22dVfZG5+5j2pinghSt+8wF0KkCVJElSo58m+0llgCpJktQxRbG5I4spjYIBqiRJUgcNo4l/XBigSpIkdUwBmw1QJUmS1CXWoEqSJKkzCuyDKkmSpG6Z3jH8BqiSJEmdU5R9UCVJktQhBZunNz41QJUkSeqaZiWp6WWAKkmS1Dlh88Irkk40A1RJkqSOKWCLTfySJEnqkmmuQd1q1BmQJEmSelmDKkmS1DHNUqfTW4NqgCpJktRBW8oAVZIkSR1hDaokSZI6pQibp3iokAGqJElSB9nEL0mSpM6wiV+SJEkdEzaXTfySJEnqiAK22AdVkiRJXWITvyRJkjqjyiZ+SZIkdcwWa1AlSZLUFc0ofmtQJUmS1Bk28UuSJKlDHMUvSZKkztnsSlKSJEnqiiJT3Qd1eksuSZKkTrIGVZIkqYO2OEhKkiRJXeE0U5IkSeqUIg6SkiRJUrc4zZQkSZI6owon6pckSVKXhC3YxC9JkqSOKKxBlSRJUsc4il+SJEmdUYQtjuKXJElSl1iDKkmSpM4oXElKkiRJnRI2O4pfkiRJXWENqiRJkjrHGlRJkiR1RlWmugZ1aCVPclKS65Jc1rNvhyRnJrmy/ffe7f4keUuSjUkuSfLwYeVLkiRpHGyurQZ+LWau+GzW8Xsl+WiSi5NcnuToFS9YH4YZmp8MHDxr3yuAz1TVnsBn2m2AQ4A929exwL8MMV+SJEnT6mTuHJ/1eiHw1araCzgA+Ick265Cvu5gaAFqVZ0N3DBr92HAKe37U4Cn9Ox/ZzXOBbZPsm5YeZMkSeqyAraQgV+LXnfu+Gz2re+ZJMA92rS3rUSZBrHafVDvX1XXAlTVtUnu1+7fGbi6J92mdt+1sy+Q5FiaWlbWcrfh5laSJGkk0leT/Rx2THJBz/b6qlo/wPlvBU4HrgHuCTyjqrYsJSPL0ZVBUnOF/DVXwvZDXg+wXXaYM40kSdI4a6aZWtIo/uurat9l3Pog4CLgscCvAWcm+UJV3byMaw5stYeHfX+m6b7997p2/yZg1550u9BE7pIkSVNpM1sN/FoBRwMfbLtdbgS+BTx4JS48iNUOUE8HjmrfHwV8pGf/c9rR/PsDP5rpCiBJkjRtirClBn+tgO8CBwIkuT/w68BVK3HhQQytiT/JqTSjv3ZMsgk4Hngd8P4kx9B8AIe3yT8OHApsBH5KE71LkiRNrS1DqEecJz7bBqCqTgD+Bjg5yaU0XTBfXlXXr3hGFjG0ALWqjpjn0IFzpC2aaQ0kSZKmXhVsXpka0VnXnTc+mzl+DfD4Fb/xgLoySEqSJEk9VqjJfiwZoEqSJHVM0wd1epc6NUCVJEnqoM19TLw/qQxQJUmSOmYZ86BOBANUSZKkzrGJX5IkSR2zxSZ+SZIkdcWwppkaFwaokiRJHWQTvyRJkjpjZqnTaTW9obkkSZI6yRpUSZKkDnKQlCRJkjrDeVAlSZLUOQ6SkiRJUnfUdA+SMkCVJEnqmMI+qJIkSeoYa1AlSZLUGQ6SkiRJUucYoEqSJKkzpn0lKQNUSZKkDnKQlCRJkrqjbOKXJElShzhISpIkSZ1jgCpJkqTOcJCUVtwZ11y8aJqDdtprFXIyHIuVb5zLBpNdPn82J7dsMNnlG+eywWSXb5LLBqMtXxmgSpIkqUumeRT/VqPOwCQa978Wl6ufmh5Jg5n23ysaX/4/YWmqHcU/6GtSGKAOwbR/Gf0fqbTypv33isaX/0/QUtjEL0mS1EH2QZUkSVKHTFaT/aAMUCVJkjrIGlRJkiR1xrSvJJWqGnUelmy77FD75cBRZ+OXljqIYRw6kE9y2WCyyzfJZQPLN5dJLhtYvi6Y5LLBypfv03Xahqrad5Br3X3PdfWQtxw9cB4uOPTvBr5XF1mDKkmS1EHTPA+qAaokSVLHFPZBlSRJUqc4il+SJEkdM8bDhJbNlaQkacK46pQ0Gaoy8GtSWIMqSRNmXEZKS5pflX1QJUmS1DH2QZUkSVKnTHMfVANUSZKkDrKJX5IkSZ1RTNagp0EZoEqSJHXQFLfwO82UJEmSumVoAWqSXZN8LskVSS5P8uJ2/w5JzkxyZfvvvdv9SfKWJBuTXJLk4cPKmyRJUqfVcOZBTXJSkuuSXLZAmgOSXNTGb2etaLn6lBrSELEk64B1VfWVJPcENgBPAZ4L3FBVr0vyCuDeVfXyJIcCfwIcCuwHvLmq9lvoHttlh9ovBw4l//1ayQmxuzh34SSXb6UnM5/k8k1y2WCyyzfJZYPJLl/XygaTXb5h/mx+uk7bUFX7DnL+2l/buXZ7/R8PfN8rD3/1gvdK8jvAj4F3VtVD5zi+PXAOcHBVfTfJ/arquoEzskxDq0Gtqmur6ivt+1uAK4CdgcOAU9pkp9AErbT731mNc4Ht2yBXI+JqNJIkjc4walCr6mzghgWSPAv4YFV9t02/6sEprFIf1CS7A/sA5wH3r6proQligfu1yXYGru45bVO7TyPStb9uJUmaJs1qUoO9VsCDgHsn+XySDUmesyJXHdDQR/EnuQfwAeAlVXVzMm90P9eBO33USY4FjgVYy91WKpuagzWokiSNRrHkeVB3THJBz/b6qlo/wPlbA78FHAjcFfhSknOr6htLycxSDTVATbINTXD67qr6YLv7+0nWVdW1bRP+TNXxJmDXntN3Aa6Zfc32Q14PTR/UoWVekiRpVApYWoB6/aD9XWfZ1F7jJ8BPkpwN7AWsaoA6zFH8AU4ErqiqN/QcOh04qn1/FPCRnv3PaUfz7w/8aKYrgEbDJn5JkkZnRE38HwEek2TrJHejGbh+xYpceQDDrEF9FPCHwKVJLmr3vQp4HfD+JMcA3wUOb499nGYE/0bgp8DRQ8yb+mATvyRJIzSEduIkpwIH0HQF2AQcD2wDUFUnVNUVST4JXAJsAd5eVfNOSTUsQwtQq+qLzN2vFJp+DbPTF/DCYeVHkiRpfAxnqdOqOqKPNP8H+D8rfvMBuJKU5mUTvyRJI1RLeE2IoY/i1/iyiV+SpBGpJY/inwgGqJrXQTvtZZAqSdKojGGNaJK39JHs5qr6i4USGKBqXgankiSN0ljWoB4GvHqRNK8ADFAlSZLGzhjWoAJvrKpTFkqQ5N6LXSS1QpNmjcJ22aH2y50mBBi61ahZHOUApUku3ySXDSa7fJNcNrB8K8GfzeGY5PKtVkvhmnVXbhh08vy77LFLrTv+Twa+13eOfsXA9+oiR/FLkiR1zcxKUoO+OiLJ/06yXZJtknwmyfVJjuz3fANUSZIkrbTHV9XNwBNplk99EPBn/Z5sH1RJkqQOGuNemNCuTkWzSuipVXVD0n8NrwGqJElSF413gPrRJF8Dfga8IMl9gVv7PdkmfkmSpC4awz6oSdYBVNUrgN8G9q2q/wJ+SjMFVV8WrEFNcvNi+QCuraoH9XtDSZIkLS7jWYN6UjuN1OeBTwJfBKiqnwA/6fciizXxf7Oq9lkoQZIL+72ZJEmS+lCMZRN/VR2SZC1wAPAHwN8n+S5NsPrJqvpuP9dZLEB9Wh/X6CeNJEmS+taNJvulqKpbaQNSgCR7AIcAb03yK1X1yMWusWCAWlVX9W4n2a73nKq6YXYaSZIkrYAxrEGdrY0dfwS8t339uJ/z+hrFn+T5wF/TjMSa+bgK+NWBcypJkqTFjXGAOl/sWFV9xY79TjP1MuA3q+r6wbMoSZKkgY1xgMoyY8d+A9Rv0kwPIEmSpGGbWep0fC0rduw3QH0lcE6S84Cfz+ysqhct9caSJEma35hOMzVjWbFjvwHq24DPApcCWwbNoSRJkgY03gHqsmLHfgPU26rquEEvLkmSpKm0rNix3wD1c0mOBT7KHatpb1jqjcfNGddcPOosaBlW8/n13uugnfZa1ftNokku32qXzZ/NlTPJZQPL1xVj3sS/rNix3wD1We2/r+zZ5zRTkiRJwzLeg6SWFTv2FaBW1R4DZkqSJElLNaZLnc5Ybuy4YICa5OFV9ZXlphmWBz3sp5xxxsWr0lS12ia9KW61y6eVM8nPbpLLNg1mnt8kPrtJ/9mc9PJNk5WKHRerQX1HkgOAheqYTwT2WeQ6kiRJGsR41qCuSOy4WIB6L2DDIjf5wSLXkCRJ0oDGdJDUisSOCwaoVbX7YHmSJEnSihjDAHWlYsd+R/FLkiRpNY1hgLpSDFAlSZI6JjW2TfwrwgBVkiSpi8Z7HtRl2aqfREk+k+TQWfvWDydLkiRJ+uVcqIO8OiTJU5O8Ick/JPmDQc7tK0AF9gBenuT4nn37DnIjSZIk9W+mmX+QV1ck+Wfgj4FLgcuA5yf5p37P77eJ/ybgQOAtST4KHDloRiVJkjSADgWcS/C7wEOrqgCSnEITrPal3xrUVNVtVfUC4APAF4H7DZpTSZIk9WEJtaddqkEFvg7s1rO9K3BJvyf3W4N6wsybqjo5yaXAC/u9iSRJkgbUrYBzUPcBrkhyfrv9COBLSU4HqKonL3RyXwFqVb1t1vYG4I8Gz6skSZL6Mt4B6quXc7LTTEmSJHVQx5rs+5LkrcB7quqs5Vyn3z6okiRJ0mKuBP4hybeTvD7J3ku5iAGqJElSF43hPKhV9eaq+m2aUfw3AO9IckWSVyd5UL/XSTv6fyztu9faOv+M3e6w76Cd9lrRe5xxzcUrer3lWsnyTXLZYLLLN8llA8u32vzZ7N8kl2+SywajLd+adVduqKqB5o9fu9Outfvzjxv4Xl//y+MGvtewJdkHOAl4WFWt6ecca1AlaRV1LQiQpGFIsk2SJyV5N/AJ4BvA0/o930FSkiRJXTSGjdxJfh84AngCcD7wXuDYqvrJINcxQJUkSeqiMQxQgVcB7wFeVlU3LPUiBqiSJEkdE8Zzmqmq+r2VuM7Q+qAmWZvk/CQXJ7k8yV+1+/dIcl6SK5O8L8m27f67tNsb2+O7DytvkiRJnTeEUfxJTkpyXZLLFkn3iCSbkzx9OUVYqmEOkvo58Niq2gvYGzg4yf7A64E3VtWewI3AMW36Y4Abq+qBwBvbdJIkSdOnmhrUQV99OBk4eKEESdbQxGFnLLscSzS0ALUaP243t2lfBTwWOK3dfwrwlPb9Ye027fEDk2RY+ZMkSeq0IdSgVtXZNPOTLuRPgA8A1y0t48s31GmmkqxJchFNAc8EvgncVFW3tUk2ATu373cGrgZoj/8IuM8w8ydJktRZI5ioP8nOwB8AJyz/aks31EFSVbUZ2DvJ9sCHgN+YK1n771y1pXf6qJMcCxwLsNvOjvGSJEmTaYmDpHZMckHP9vqqWj/A+W8CXl5Vm0fZkL0qEV5V3ZTk88D+wPZJtm5rSXcBrmmTbQJ2BTYl2Rq4F3NUQbcf8npoVpJahexLkiStvqVFOdcvcyWpfYH3tsHpjsChSW6rqg8v45oDG+Yo/vu2NackuSvwOOAK4HPAzIiwo4CPtO9Pb7dpj3+2xnkdVkmSpKVaSvP+CkRNVbVHVe1eVbvTjAl6wWoHpzDcGtR1wCntSLCtgPdX1ceSfJUmMn8NcCFwYpv+ROBdSTbS1Jw+c4h5kyRJ6rRhzIOa5FTgAJquAJuA42kGslNVI+132mtoAWpVXQLsM8f+q4BHzrH/VuDwYeVHkiRprAwhQK2qIwZI+9yVz0F/Ms6t6PvutbbOP2O3eY8ftNNeS7ruGddcvNQsrZpJLhtYvvmMQ/kmuWww2eWb5LKB5ZvLJJcNulO+Neuu3DBov9C7/squ9cAjjxv4Xpf9w3ED36uLHAYvSZLUReNbh7hsBqiSJElds0KDnsaVAaokSVLHhLkniJ8WBqiSJEldNMU1qENd6lTD05WO35Kk8bXUwUfSsFmDOqb8pSJJWi4rO7ptGPOgjgsDVEmSpC6a4gDVJv4x5V+9kiRNuBEsddoV1qCOKZv4JUmaYGUTvyRJkrrGAFWSJEldYg2qJEmSusUAVZIkSV0yzTWoqRrf0u+719o6/4zd+kq72KCicR4V38+AqUku3ySXDSa7fONcNpjs8k1y2WCyy+fvle6Vbc26KzdU1b6DnHO3++5aD37acQPf68K3HTfwvbrIGlRJkqQuGt86xGVzHlRJ0sCc6k4artA08Q/6mhTWoEqSBtbFZlRp4kxQwDkoA1RJkqQOyhiPE1ouA1RJkqSumbClSwdlgCpJktRBk9SndFAGqJIkSV00xQHq1MyDOqN35OkkdvKflvJNctnA8o2bSS4bTE/5JrlsMNnl63rZljIP6t133LUe8qQ/HfheF5z8UudBlSRJ0nBMcxP/1NWgSpIkraal1qD+5hMGr0H98jutQZUkSdIwTNjE+4MyQJUkSeoiA1RJkiR1xcxSp9PKAFWSJKmLxnic0HIZoEqSJHWQNaiSJEnqDpc6lSRJUtdky6hzMDoGqJIkSV1kDaokSZK6xD6okiRJ6o7CUfySJEnqFmtQJUmS1C0GqJIkSeoKV5KSJElSt1RNdR/UrUadAUmSJKmXNaiSJEkdZBO/JEmSusUAVZIkSV1iDaokSZK6o4At0xuhGqBKkiR10fTGpwaokiRJXTTNTfxDn2YqyZokFyb5WLu9R5LzklyZ5H1Jtm3336Xd3tge333YeZMkSeqsmblQB3ktIslJSa5Lctk8x5+d5JL2dU6SvVa8XH1YjXlQXwxc0bP9euCNVbUncCNwTLv/GODGqnog8MY2nSRJ0lRKDf7qw8nAwQsc/xbwu1X1MOBvgPXLLsgSDDVATbIL8ATg7e12gMcCp7VJTgGe0r4/rN2mPX5gm16SJGm61BJfi1226mzghgWOn1NVN7ab5wK7LLUIyzHsGtQ3AX8ObGm37wPcVFW3tdubgJ3b9zsDVwO0x3/Upr+DJMcmuSDJBT/44eZh5l2SJGkkAqRq4Bew40yc1L6OXUY2jgE+sSIFGtDQBkkleSJwXVVtSHLAzO45klYfx27fUbWetrp5373WTnH3YUmSNNG2LJ5kDtdX1b7LvXWS36MJUB+93GstxTBH8T8KeHKSQ4G1wHY0NarbJ9m6rSXdBbimTb8J2BXYlGRr4F4sUAUtSZI0ydLHoKeh3Dd5GE33zEOq6oejyMPQmvir6pVVtUtV7Q48E/hsVT0b+Bzw9DbZUcBH2vent9u0xz9bNaInI0mSNEpD6oO6mCS7AR8E/rCqvrH8Ky7NKOZBfcfb7AIAAAkxSURBVDnw3iSvAS4ETmz3nwi8K8lGmprTZ44gb5IkSR3Q37RRg0pyKnAATV/VTcDxwDYAVXUC8GqaMUD/3I5Vv20lugwMalUC1Kr6PPD59v1VwCPnSHMrcPhq5EeSJKnrhjFRf1Udscjx5wHPW/k7D8aVpCRJkrpoins6rsZE/ZIkSVLfrEGVJEnqmoIsbZqpiWCAKkmS1EVT3MRvgCpJktRF0xufGqBKkiR10agm6u8CA1RJkqQuMkCVJElSZxTgIClJkiR1RSib+CVJktQxBqiSJEnqFANUSZIkdYZ9UCVJktQ19kGVJElStxigSpIkqTvKAFWSJEkdUhigSpIkqWMcJCVJkqQumeZBUluNOgOSJElSL2tQJUmSumiKa1ANUCVJkrqmgC0GqJIkSeoMp5mSJElS1xigSpIkqVMMUCVJktQZ9kGVJElStxTU9M7Ub4AqSZLURTbxS5IkqTNs4pckSVLnWIMqSZKkTjFAlSRJUnc4Ub8kSZK6pIAtjuKXJElSl1iDKkmSpE4xQJUkSVJ3lNNMSZIkqUMKaopXktpq1BmQJEmSelmDKkmS1EU28UuSJKlTHCQlSZKkzqhyHlRJkiR1jDWokiRJ6pKyBlWSJEndUdagSpIkqUMKR/FLkiSpY5yofziSfDvJpUkuSnJBu2+HJGcmubL9997t/iR5S5KNSS5J8vBh5k2SJKmrCqgtNfBrMUlOSnJdksvmOd6JeGw1VpL6varau6r2bbdfAXymqvYEPtNuAxwC7Nm+jgX+ZRXyJkmS1D1VTQ3qoK/FnQwcvMDxTsRjo1jq9DDglPb9KcBTeva/sxrnAtsnWTeC/EmSJI3cMGpQq+ps4IYFknQiHht2gFrAp5JsSHJsu+/+VXUtQPvv/dr9OwNX95y7qd0nSZI0fYZTg7qYTsRjwx4k9aiquibJ/YAzk3xtgbSZY9+d/hRoA92ZYPfna9ZdOWcfCnXejsD1o86ElsRnN958fuPLZze+fn3QE27hxjM+XaftuIR7rZ0Z99NaX1XrBzi/r3hs2IYaoFbVNe2/1yX5EPBI4PtJ1lXVtW2V8XVt8k3Arj2n7wJcM8c11wPrAZJc0NO3VWPEZze+fHbjzec3vnx242tWwNiXqlqon+gw9RWPDdvQmviT3D3JPWfeA48HLgNOB45qkx0FfKR9fzrwnHb02P7Aj2a6AkiSJGlVdCIeG2YN6v2BDyWZuc97quqTSb4MvD/JMcB3gcPb9B8HDgU2Aj8Fjh5i3iRJkqZOklOBA4Adk2wCjge2AaiqE+hIPDa0ALWqrgL2mmP/D4ED59hfwAsHvM0gfSrULT678eWzG28+v/HlsxtfnXl2VXXEIseXEo+tuNQUr/MqSZKk7hnFPKiSJEnSvMY2QE1ycJKvt0txvWLxM7Sakuya5HNJrkhyeZIXt/td6nYMJFmT5MIkH2u390hyXvvc3pdk23b/Xdrtje3x3UeZb0GS7ZOcluRr7ffvt/3ejYckf9r+vrwsyalJ1vrd66a5lgtdyvcsyVFt+iuTHDXXvabVWAaoSdYA/0SzHNdDgCOSPGS0udIstwEvrarfAPYHXtg+I5e6HQ8vBq7o2X498Mb2ud0IHNPuPwa4saoeCLyxTafRejPwyap6MM04gCvwe9d5SXYGXgTsW1UPBdYAz8TvXledzJ2XCx3oe5ZkB5oBSvvRTMN5/ExQqzENUGke5MaquqqqfgG8l2ZpLnVEVV1bVV9p399C8z/JnXGp285LsgvwBODt7XaAxwKntUlmP7eZ53kacGCbXiOQZDvgd4ATAarqF1V1E37vxsXWwF2TbA3cDbgWv3udNM9yoYN+zw4CzqyqG6rqRuBM7hz0Tq1xDVA7sQyX+tM2Pe0DnIdL3Y6DNwF/DsysmXcf4Kaquq3d7n02v3xu7fEftek1Gr8K/AB4R9tF4+3tPNR+7zquqr4H/D3N9IvX0nyXNuB3b5wM+j3z+7eAcQ1QO7EMlxaX5B7AB4CXVNXNCyWdY5/PdJUleSJwXVVt6N09R9Lq45hW39bAw4F/qap9gJ9wezPjXHx+HdE27R4G7AHsBNydpml4Nr9742e+Z+UzXMC4BqidWIZLC0uyDU1w+u6q+mC7+/szTYhZwlK3GrpHAU9O8m2arjOPpalR3b5tdoQ7PptfPrf2+L24c7OXVs8mYFNVnddun0YTsPq9677HAd+qqh9U1X8BHwT+G373xsmg3zO/fwsY1wD1y8Ce7ejGbWk6kp8+4jypR9sX6kTgiqp6Q88hl7rtsKp6ZVXtUlW703yvPltVzwY+Bzy9TTb7uc08z6e36a0BGJGq+k/g6iS/3u46EPgqfu/GwXeB/ZPcrf39OfPs/O6Nj0G/Z2cAj09y77YG/fHtPjHGE/UnOZSmZmcNcFJVvXbEWVKPJI8GvgBcyu19GV9F0w/1/cButEvdVtUN7S/kt9J0EP8pcHRVXbDqGdcvJTkAeFlVPTHJr9LUqO4AXAgcWVU/T7IWeBdNH+MbgGe2q8hpRJLsTTPAbVvgKpplCrfC713nJfkr4Bk0s6BcCDyPpk+i372OSc9yocD3aUbjf5gBv2dJ/ojm/40Ar62qd6xmObpsbANUSZIkTaZxbeKXJEnShDJAlSRJUqcYoEqSJKlTDFAlSZLUKQaokiRJ6hQDVEmSJHWKAaqkqZZk9yQ/S3LRgOc9I8nGJB8bVt4kaVoZoEoSfLOq9h7khKp6H81E6pKkFWaAKmliJXlEkkuSrE1y9ySXJ3noIufsnuRrSd6e5LIk707yuCT/keTKJI9crfxL0rTaetQZkKRhqaovJzkdeA1wV+DfquqyPk59IHA4cCzwZeBZwKOBJ9MsS/iU4eRYkgQGqJIm31/TBJm3Ai/q85xvVdWlAEkuBz5TVZXkUmD3oeRSkvRLNvFLmnQ7APcA7gms7fOcn/e839KzvQX/sJekoTNAlTTp1gP/C3g38PoR50WS1AdrAiRNrCTPAW6rqvckWQOck+SxVfXZUedNkjS/VNWo8yBJI5Nkd+BjVbXg6P55zj0AeFlVPXGFsyVJU80mfknTbjNwr6VM1A/8M3DjUHIlSVPMGlRJkiR1ijWokiRJ6hQDVEmSJHWKAaokSZI6xQBVkiRJnWKAKkmSpE75vz32CisKU0UFAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Create model with finer discretization\n",
    "\n",
    "# Define new FD grid parameters\n",
    "dh1 = 5.0\n",
    "Xmax = 1000.0\n",
    "Zmax = 500.0\n",
    "\n",
    "# Calculate coordinates of grid points\n",
    "NX1, NZ1, x1, z1, XX1, ZZ1 = coord_def(Xmax,Zmax,dh1)\n",
    "\n",
    "# Define model with finer discretization\n",
    "vp_sine1 = create_sin_model(NX1,NZ1,dh1)\n",
    "\n",
    "# Plot Vp-model\n",
    "cax = plt.pcolormesh(x1, z1, vp_sine1)\n",
    "plt.title(\"Sine layer P-wave velocity model\" )\n",
    "plt.xlabel(\"x [m]\")\n",
    "plt.ylabel(\"z [m]\")\n",
    "cbar = plt.colorbar(cax, orientation='vertical', pad=0.02)\n",
    "cbar.set_label('Vp [m/s]', labelpad=10)\n",
    "plt.gca().invert_yaxis()\n",
    "# plt.savefig('vp_sine1.pdf', bbox_inches='tight', format='pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks better, but the interface is still quite blocky. This is a serious disadvantage of  the Cartesian grid and therefore also for the finite-difference approach. We get these \"stair-case\" discretization errors, which become a problem in seismic modelling, because they lead to artificial diffractions. Try to minimize the stair-case artifacts by using an even finer model discretization."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## We learned:\n",
    "\n",
    "* Creating 2D Cartesian grids and P-wave velocity models\n",
    "* Visualization of discrete models\n",
    "* The Cartesian grid has serious limitations to accurately discretize complex layer interfaces which not coincide with the Cartesian coordinate axes"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
